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I.  INTRODUCTION 


This  report  discusses  the  work  accomplished  during  the  Air  Force  Office  of 
Scientific  Research  Grant  AFOSR-86-0142.  There  were  three  primary  objectives  of  this 
work.  The  first  was  to  develop  a  method  for  estimating  the  stress  intensity  factors  induced 
by  temperature  gradients  in  thermo-mechanical  fatigue.  The  second  was  to  use  this  method 
to  estimate  the  importance  of  the  thermal  contribution  relative  to  the  mechanical 
contribution.  To  do  this,  we  take  into  account  several  factors,  for  example,  the  dependence 
of  the  stress  intensity  on  the  frequency,  the  crack  geometry  and  the  mode  of  thermal 
loading  (i.e.  heat  lamp  vs.  electric  resistance  heating).  The  final  goal  was  to  extend  the 
developed  method  to  composite  material  applications  and  perform  a  preliminary  study  to 
assess  the  relative  importance  of  thermal  transient  effects  in  composite  materials.  All  these 
objectives  have  been  met  and  a  detailed  discussion  follows. 

Gas  turbine  engine  components  typically  undergo  complex  mechanical  and  thermal 
operating  conditions  which  may  give  rise  to  crack  propagation  and  eventual  fatigue  failure. 
Although  the  role  of  mechanical  loading  in  fatigue  has  been  investigated,  it  is  of 
fundamental  imponance  to  understand  how  cyclic  thermal  conditions  affect  crack  growth 
and  fatigue  failure.  This  study  is  an  attempt  to  estimate  the  importance  of  cyclic  thermal 
loads  in  theimo-mechanical  fatigue. 

Laboratory  tests  are  currently  used  to  assess  the  safe  thermomechanical  fatigue 
limits  of  engine  components.  However,  there  are  three  aspects  of  the  thermal  loading 
which  potentially  vary  widely  between  laboratory  simulations  and  actual  in-service 
conditions.  These  are  the  frequency  of  the  thermal  cycle,  the  symmetry  of  the  temperature 
field  about  the  crack  axis,  and  the  geometry  of  the  initial  crack. 

The  frequency  of  the  thermal  cycle  is  related  to  the  heating  rate  which  may  be 
significantly  higher  or  lower  in  engine  environments  than  in  fatigue  testing  situations. 
Currently,  heating  rates  in  fatigue  testing  are  limited  to  a  level  which  is  believed  to  not 
induce  significant  thermal  stresses  at  the  crack  tip.  Consequently,  it  is  unknown  how 
accurately  results  of  stress  intensity  testing  can  be  applied  to  engine  fatigue  standards. 
Determining  the  variation  of  the  stress  field  with  frequency  provides  a  means  of  accounting 
for  the  high  heating  rates  in  both  engine  components  and  lab  specimens  thus  improving  the 
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validity  and  accuracy  of  the  test  method.  Further,  having  a  means  to  account  for  rapid 
thermal  cycling  also  may  provide  an  economic  benefit  in  that  it  may  be  possible  to  increase 
the  rate  at  which  specimens  are  cycled. 

Test  specimens  are  usually  symmetrically  heated  with  respect  to  the  crack  axis. 
This  is  a  fundamental  difference  from  what  happens  to  engine  components  in  that  cracks 
are  usually  subjected  to  temperature  gradients  that  may  induce  significant  thermal  stresses  at 
the  crack  tip.  The  reason  is  related  to  the  thermal  boundary  conditions  specified  along  the 
axis  of  the  crack  tip.  The  symmetry  condition  on  this  axis  is  identical  to  the  insulation 
condition  of  the  open  crack;  the  normal  gradient  of  the  temperature  is  zero.  The 
antisymmetry  condition,  on  the  other  hand,  is  that  the  temperature  itself  is  zero  along  this 
boundary.  This  gives  rise  to  a  discontinuity  in  the  boundary  condition  at  the  crack  tip  and 
results  in  a  thermal  field  which  changes  as  the  crack  grows.  The  effect  of  this  thermal 
field/crack  interaction  is  currently  unaccounted  for  in  fatigue  analysis. 

The  effect  of  geometry  on  stress  intensity  is  also  being  investigated.  Cracks  can  be 
located  on  the  specimen  edge  or  in  its  center.  The  important  distinction  in  the  crack  location 
is  that  the  edge  cracked  component  corresponds  to  a  simply-connected  region  while  the 
center  cracked  component  corresponds  to  a  multiply-connected  region.  Thus,  an  additional 
aspect  to  be  investigated  is  the  effect  of  geometry  on  stress  intensity. 

The  means  for  evaluating  the  importance  of  each  of  these  phenomena  is  the  stress 
intensity  factor,  K.  The  stress  intensity  factor  is  a  quantitative  measure  of  the  strength  of 
the  elastic  stress  field  near  the  crack  tip  and  is  of  interest  because  of  the  general  belief  that  it 
is  a  dominant  parameter  in  controlling  crack  instability.  We  can  use  the  standard  definitions 
of  stress  intensity  factor  for  thermoelastic  problems  since  it  has  been  shown  that  the  local 
behavior  of  the  thermal  stresses  near  the  crack  tip  is  the  same  as  the  behavior  of  mechanical 
stresses,  that  of  the  singularity  [1]. 

The  first  work  dealing  with  two-dimensional  stress  fields  associated  with  cracks 
was  published  by  Florence  and  Goodier  [2]  and  Sih  [1].  The  method  of  complex 
representation  of  the  elastic  state  was  used  to  derive  the  theoretical  stress  intensity  factors  in 
infinite  two-dimensional  cracked  regions  subjected  to  uniform  heat  flow  for  both  the 
symmetric  and  antisymmetric  problems.  The  expressions  derived  in  [1]  for  the 
thermoelastic  crack-tip  stress  field  are  identical  to  those  given  by  Irwin  [3]  for  traction 
boundary  condition  problems.  It  is  noted  that  the  real  parts  of  the  crack-tip  stress  field  are 
identical  to  the  singular  terms  of  the  Williams  [4]  eigenfunction  relations. 

Several  investigators  have  used  stress  intensity  factors  in  studying  problems  of 
thermoelastic  stress  intensity  in  cracked  regions  subjected  to  nonuniform  temperature 
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fields.  Early  work  of  Konishi  and  Atsumi  [5]  and  of  Sekine  [6,7,8]  used  singular  integral 
equations  to  estimate  the  normal  thermal  stress  components  on  the  crack  plane  in  semi- 
infinite  plates  with  uniform  heat  flow  perpendicular  to  the  crack.  From  these  solutions  the 
stress  intensity  factors  were  evaluated.  This  method  has  also  been  used  by  Shall  [9]  and 
Das  [10]  for  determining  stress  intensity  factors  in  an  externally  cracked  infinite  solid  for 
both  even  and  odd  temperature  distributions,  by  Das  [11]  for  cracked  infinite  cylinders  with 
constant  temperature  and  by  Herrmann  and  Kuemmerling[12]  for  cracked  infinite  cylinders 
with  a  point  source  of  heat. 

Because  of  the  finite  domains  of  the  present  analysis,  however,  solution  by  these 
methods  is  not  appropriate  and  an  alternative  way  to  evaluate  the  stress  intensity  factor  from 
numerically  defined  stress  fields  is  required.  Path  independent  integrals  have  been 
developed  to  estimate  the  stress  intensity  factors  at  cracks  and  notches  when  numerical 
approaches,  such  as  the  finite  element  method,  are  used  to  compute  singular  stress 
distributions.  The  path  independent  integrals  offer  the  distinct  advantage  that  the  errors  in 
the  analysis  stem  only  from  the  numerical  approximations  used  and  not  from  additional 
sources  [13]  and  so  offer  a  relatively  accurate  result  with  modest  computational  effort. 

Several  path  independent  integrals  for  thermoelastic  stress  problems  have  been 
developed.  The  integral  of  Gurtin  [14]  consists  of  the  well-known  J-integral  for  isothermal 
stress  fields  and  three  additional  terms  which  relate  to  the  temperature  field.  However,  the 
development  of  the  integral  relies  on  specifying  the  temperature  field  to  be  symmetric  about 
the  crack  axis  and  zero  along  the  crack  surfaces,  requirements  that  are  not  in  general  met. 
Wilson  and  Yu  [15]  have  developed  an  integral  similar  to  the  J-iniegral  which  consists  of 
an  additional  area  integral  containing,  not  only  the  temperature  field  but  also  the  gradient  of 
the  temperature  field.  Aoki  et  al.  [16,17,18]  have  derived  another  set  of  path  independent 
integrals  for  thermoelastic  problems.  As  in  the  integral  of  Wilson  and  Yu,  each  of  these 
integrals  require  area  integration  of  the  gradient  of  the  temperature  distribution,  implying 
that  a  highly  accurate  temperature  solution  near  the  crack  tip  is  necessary  in  order  to 
evaluate  the  stress  intensity  factor  precisely.  Finally,  Kuo  and  Riccardella  [19]  have 
developed  line  integrals  for  thermoelastic  stress  fields  which  eliminate  the  need  for  area 
integration  but  in  doing  so,  restrict  the  problem  to  one  where  no  body  forces  or  distributed 
heat  sources  are  permitted. 

The  appropriate  integral  for  our  work  follows  from  the  development  of  the  H- 
integral  by  Sinclair,  Okajima  and  Griffin  [13]  and  takes  into  account  the  Duhamel- 
Neumann  analogy  (for  example,  [20])  for  modeling  the  thermoelastic  problem.  The  path- 
independent  integral  consists  of  the  H-integral  and  an  additional  area  integral  containing  the 
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temperature  field.  Although  the  derived  integral  loses  some  computational  efficiency 
through  calculation  of  the  area  integral,  it  does  not  require  the  gradient  of  the  temperature 
distribution.  Since  our  temperature  distributions  are,  in  general,  numerically  determined, 
calculating  thermal  gradients  would  be  prohibitively  inaccurate. 

In  our  approach  to  studying  thermal  fatigue  loading,  we  concentrate  on  the  effect 
that  the  thermal  cycle  has  on  thermally  induced  crack  tip  stresses.  Since  any  periodic 
function  can  be  represented  exactly  as  a  Fourier  series  of  harmonic  functions,  we  assume 
all  thermal  cycles  are  harmonic  in  time.  By  this  assumption,  we  can  consider  the  derived 
solutions  to  the  harmonic  thermal  cycles  as  Fourier  components  of  the  true  solution  to 
general  thermal  cycling  cases.  The  actual  solution  may  be  approximated  by  its  fundamental 
harmonic  or  as  closely  as  desired  by  adding  higher  frequency  components. 

Considering  only  sinusoidal  thermal  loadings  provides  the  advantage  that  we  can 
eliminate  time  as  an  explicit  variable.  Without  time  as  a  variable,  there  is  no  need  to  track 
stress  fields  as  a  function  of  time.  Instead,  its  role  is  replaced  by  the  parametrical 
dependence  of  frequency.  With  this,  we  can  calculate  stress  intensity  factors  for  sinusoidal 
thermal  loadings  as  a  function  of  frequency  of  the  thermal  cycle. 

The  thermal  loadings  we  consider  are  selected  after  careful  evaluation  of  fatigue 
testing.  We  choose  simple  loading  conditions  which,  when  combined  through  linear 
superposition,  describe  more  complex  loading  cases.  The  goal  in  our  approach  is  to 
develop  stress  intensity  solutions  to  simple  problems  that  can  be  combined  later  to  represent 
stress  intensities  of  a  large  variety  of  fatigue  loading  cases.  With  this,  we  have  outlined  an 
approach  to  calculating  thermal  stress  intensities  as  functions  of  frequency  in  fatigue 
loading.  We  evaluate  the  importance  of  the  stress  intensity  by  comparing  it  to  those 
induced  in  mechanical  fatigue,  for  example,  the  results  of  Wilson  and  Warren  [21]. 

We  begin  in  Section  2  by  addressing  the  basic  problem  of  a  symmetrically  and 
uniformly  heated  center-cracked  fatigue  specimen,  presenting  the  formulation,  a  discussion 
of  the  solution  methods,  and  the  results  of  this  simple  analysis.  In  Section  3,  we  consider 
extensions  of  the  basic  problem.  First,  we  investigate  some  geometrical  variations, 
considering  both  edge-cracked  regions  and  crack-  and  sample-length  scaling.  Secondly, 
we  address  the  issue  of  antisymmetric  fatigue  heating  and  the  importance  of  Mode  II  stress 
intensity  factors.  Because  specimen  heating  is  often  accomplished  by  applying  an  electric 
potential  to  the  specimen,  we  investigate  this  mode  of  heating  in  Section  4.  Evaluating 
electric  potential  heating  is  important  because  it  induces  a  singular  heat  source  at  the  crack 
tip  which  may  serve  to  give  a  much  different  dependence  on  frequency  than  the  regular 
heating  fields  studied  in  Sections  2  and  3.  In  Section  5,  we  close  the  analysis  of  the 
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homogeneous  linear  elastic  material  by  investigating  the  implications  the  results  obtained 
have  in  actual  fatigue  testing. 

Finally,  in  the  interest  of  understanding  the  effects  of  thermal  transients  in 
composite  materials,  we  extend  our  approach  of  fatigue  analysis  to  bimaterial  composites. 
In  Section  6  we  discuss  the  fundamentals  of  a  pilot  study  performed  to  investigate  the 
rudimentary  effects  of  diffusivity  and  thermal  expansion  mismatches  in  bimaterial 
composites  in  which  cracks  occur  normal  to  the  interface.  This  particular  geometry  was 
selected  for  initial  investigation  because  of  its  physical  interest:  it  can  correspond  to  the 
problem  of  a  crack  turning  at  or  penetrating  the  interface. 

The  intent  of  the  composite  analysis  is  to  perform  a  parametrical  study  to  establish 
the  importance  of  differences  in  thermal  diffusivity  relative  to  differences  in  the  thermal 
expansion  and  elastic  modulus  for  composite  materials.  These  effects  may  be  especially 
important  in  materials  with  low  diffusivity  such  as  ceramic  comjjosites.  Our  resulting 
improvement  in  understanding  these  mechanisms  should  allow  the  development  of 
composites  with  improved  durability  for  high  temperature  applications. 


n.  TRANSIENT  UNIFORM  HEAT  SOURCE  IN  A  CENTER  CRACKED  PLATE 
Formulation 

In  order  to  investigate  the  effects  of  transient  thermal  loadings  on  crack  tip  stress 
fields,  we  formulate  the  following  problem.  Our  analysis  is  directly  applicable  to  fatigue 
testing  so  we  choose  a  geometry  which  represents  typical  fatigue  test  specimens.  This  is 
the  center-cracked  plate  pictured  in  Figure  1.  Geometrical  symmetry  permits  analysis  of 
only  the  upper  right  quadrant  of  the  center-cracked  plate,  the  region  shaded  in  the  figure. 
For  simplicity,  we  assume  the  plate  to  be  in  a  state  of  plane  strain.  The  region,  R  ,  of  our 
analysis,  then,  is  rectangular  in  the  x,  y  plane  and  is  deHned  by : 

R  =  (x,y  I  O^xSW,  O^y^L) 


where  R  is  assumed  to  be  homogeneous,  isotropic  and  linear  elastic.  The  boundary  of  R 
is  denoted  by  9/?  and  contains  the  crack  of  length  o  along  the  line 

0<x<a,y  =  0 

We  seek  then,  the  two-dimensional  thermoelasdc  stress  vector  a  =  (Ox  ,Oy  ,txy ) 
and  displacement  vector  u  =  (u,  v)  satisfying  the  plane-stra.'n  uncoupled  quasi-static 
equations  of  thermoelasticity.  These  equations  are  derived  by  neglecting  the  effects  of 
elasticity  on  the  temperature  and  ignoring  the  inertial  term  in  the  equation  of  motion.  The 
result  is  the  degeneration  of  heat  conduction  and  elasticity  into  two  separate  problems. 
With  these  simplifying  assumptions,  the  temperature  field  is  determined  solely  by 
conduction  and  the  stress  and  displacement  fields  are  computed  for  each  instantaneous 
temperature  distribution  according  to  the  equations  of  linear  thermoelasticity. 

The  temperature  distribution  ,T,  varies  with  time,  t,  and  satisfies  the  planar  heat 
conduction  equation  in  R  : 

T  (x,y,t )  =  -^  T«  (x,y,t )  -  H  ( x,y,t ) 
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where  is  the  Laplacian  operator  in  two  dimensions,  D,  the  thermal  diffusivity,  is  equal 
to  the  thermal  conductivity  divided  by  the  heat  capacity  and  the  density,  and  the  subscript 
following  a  comma  denotes  differentiation  with  respect  to  that  variable.  The  forcing 
function  H(x,y,t)  is  a  net  distributed  heat  source  per  unit  area  term  which  takes  into  account 
the  effects  of  radiation  and  convection  from  the  surface.  The  boundary  conditions  imposed 
on  the  temperature  field  are  that  the  plate  is  insulated  along  the  right  edge  and  along  the 
crack  faces,  and  the  temperature  is  symmetric  along  the  y-axis  and  along  the  uncracked 
portion  of  the  x-axis.  Since  the  conditions  of  insulation  and  symmetry  are  identical  (the 
normal  temperature  gradient  equals  zero)  these  conditions  may  be  stated: 


T,j  =0  on  X  =  W 

(0<y<L) 

T,j  =0  on  X  =  0 

(0<y<L) 

(4) 

T,y  =  0  on  y  =  0 

(0<x<W) 

Along  the  remaining  edge  of  the  boundary,  a  spatially  constant  temperature  is  imposed: 

T(x,y,t)  =0o(t)  ony  =  L(0<x<W)  ^5^ 

The  equations  of  uncoupled  quasi-static  linear  thermoelasticity  are  identical  to  those 
of  linear  elasticity  with  the  exception  of  the  stress  displacement  law.  This  being  true,  the 
stress  field  obeys  the  equations  of  equilibrium  in  the  absence  of  a  body  force  field,  which 
hold  in  the  region  R  : 

CJjj  ,  ,  X^y  ,  y  =  0 

S*y  =  0  (7) 

The  stress  and  displacement  vectors  must  also  satisfy  the  plane-strain  stress-displacement 
relations  in  R : 

=  0~-^v)  ^(l-'’)u,x  +  vv,y]  -  PT 
2H 

"  TT-iv)  [(i-v)v.y  +  V  u.,  ]  -  PT 
^xy  '^»x) 


(8) 
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where  P  combines  material  constants  in  the  form,  P  =  a  E  /  (1  -  2v  )  for  the  plane  strain 
case  (  P  =  a  E  /  (1  -  V  )  for  the  plane  stress  case  )  and  the  constants  a,E,v,  and  p. 
represent  the  coefficient  of  linear  thermal  expansion,  the  elastic  modulus,  Poisson's  ratio, 
and  the  shear  modulus,  respectively. 

The  boundary  conditions  for  the  thermoelastic  problem  are;  traction-free  conditions 
on  the  upper  and  right  edges  of  the  region  and  along  the  crack  face, 

Oy  =  0  ,  Txy  =  0  on  y  =  L  ( 0  <  X  <  W) 

=  0 ,  Txy  =  0  on  X  =  W  ( 0  <  y  <  L) 

Oy  =  0,Xxy  =  0  ony  =  0(0<x<a) 

and  symmetry  conditions  along  the  y-axis  and  along  the  x-axis  in  the  region  not  containing 
the  crack. 


u  =  0  and  Xj^y  =0  on  x  =  0  (0<y<L) 

V  =  0  and  x^y  =  0  on  y  =  0  ( a  <  x  <  W ). 


(10) 

(11) 


The  means  for  evaluating  the  importance  of  transient  temperature  fields  on  crack  tip 
stress  fields  is  the  stress  intensity  factor,  K.  The  stress  intensity  factor  is  a  quantitative 
measure  of  the  strength  of  the  elastic  stress  field  near  the  crack  tip  and  is  defined  for  the 
crack-opening  mode.  Mode  I,  by  the  relation 

K  =  lim  -J  27t(x-a )  Oy  on  y  =  0 

(12) 


and  for  the  edge-sliding  mode.  Mode  II,  by  the  relation 
K=  lim  J2n(x-a )  x,v  on  y  =  0 


(13) 


Considering  transient  temperature  fields  typically  implies  that  the  stress  and 
displacement  vectors,  and  consequently,  the  stress  intensity  factor  must  be  computed  at 
instantaneous  points.  Herein  we  impose  a  fundamental  condition;  we  consider  only 
temperature  fields  which  vary  sinusoidally  in  time.  In  essence,  the  temperature  field, 
T(x,y,t),  is  taken  as  the  sum  of  harmonic  components  in  which  O)  is  the  frequency  of  the 
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excitation  and  Tg  and  Tc  are  the  spatially  varying  magnitudes  of  their  respective 
components: 


T  ( x,y,t )  =  Tj  (x,y)  sin  cot  +  (x,y)  cos  cot 


(14) 


By  representing  the  temperature  in  this  form  and  considering  only  harmonic  heat  sources, 
H(x,y,t)  =  h(x,y)  cos  cot  ^  we  benefit  in  two  ways.  First,  the  need  for  instantaneous 
computations  is  eliminated  since  time  is  no  longer  an  explicit  variable.  This  is  because  we 
calculate  the  stress  and  displacement  fields  induced  by  each  of  the  spatial  functions,  Ts(x,y) 
and  Tc(x,y),  independently.  Then,  if  Kg  and  are  the  resultant  stress  intensity  factors 
asscx:iated  with  these  fields,  respectively,  the  total  stress  intensity  factor  induced  by  the 
temperature  field  T(x,y,t)  is  also  obtained  by  hanoonic  sum: 

K(co,t)  =  K,(co)  sin  cot  +  Kc(co)  coscot 


Secondly,  because  any  pericxiic  function  can  be  represented  exactly  by  Fourier 
series  of  sinusoidal  functions,  solutions  derived  with  this  approach  may  be  considered  as 
Fourier  components  of  a  solution  to  a  more  complex  periodic  temperanue  profile.  This  is 
imponant  in  fatigue  studies  where  actual  thermal  cycling  experiments  may  not  always 
utilize  sinusoidal  temperature  fields  and  the  in-service  conditions  the  experiments  are 
intended  to  simulate  may  also  be  more  complex  than  harmonic  functions. 

The  stress  intensity  factors  in  (15)  are  denoted  as  functions  of  the  frequency. 
Eliminating  time  explicitly  requires  its  role  to  be  replaced  by  the  parametrical  dependence  on 
©.  This  is  because  substitution  of  (14)  into  (3)  results  in  coupled  equations  for  Tg  and  Tc 
in  which  frequency  appears  as  a  parameter 


V^T,  =  -gT. -h{x,y) 


V^T, 


T 


(16) 


A  range  of  frequencies  must  therefore  be  analyzed  before  a  general  solution  can  be 
obtained. 


^If  the  heat  source  has  a  sine  component  as  well,  its  contribution  may  be  readily  computed  by  introducing  a 
90  degree  phase  shift  in  the  solution  computed  here. 
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Problem  definition 

The  goal  of  this  analysis  is  to  form  a  basis  for  understanding  how  transient 
temperature  fields  affect  stress  fields  in  a  region  containing  a  crack.  A  desired  result  would 
be  to  apply  solutions  to  fatigue  testing  cases  so  that  we  can  make  estimates  of  the 
importance  of  the  thermal  field  with  respect  to  mechanically  induced  stresses  the  specimen 
may  experience.  The  trouble  with  this  is  that  many  different  variations  in  loadings  are 
generally  used  and  it  would  be  impossible  to  explicitly  calculate  temperature  fields 
corresponding  to  each  case.  Instead,  we  choose  to  use  the  linearity  of  the  stress  and 
displacement  equations  and  hence  their  vector  solutions  in  defming  the  loading  conditions 
to  be  investigated.  By  carefully  choosing  speciHc  thermal  loading  cases  and  obtaining  their 
stress  field  solutions,  superposition  may  be  used  to  represent  a  variety  of  cases  by  linear 
combination  of  the  conditions  and  so  more  general  solutions  may  be  obtained. 

In  order  to  define  the  appropriate  loading  conditions,  we  examine  fatigue  testing  a 
little  more  closely.  In  general,  two  methods  of  heating  are  utilized  to  cycle  specimen 
temperatures.  These  are  a  heat  lamp  with  forced  air  cooling,  and  direct  resistance  heating 
with  grip  cooling.  Grips  hold  the  specimens  in  place  and  often  extend  the  width  of  the 
specimen  at  distances  considered  remote  from  the  crack.  In  order  to  account  for  heat 
conduction  out  of  the  crack  region,  two  alternatives  are  considered.  First,  it  would  be 
possible  to  consider  conduction  conditions  on  these  ends  in  which  the  heat  flux  is 
proportional  to  the  temperature  at  the  ends.  However,  this  approach  would  require  varying 
the  proportionality  constant  depending  on  the  grip  configuration  and  would  consequently 
require  a  number  of  simulations  in  order  to  represent  the  entire  range  of  possible 
conduction  conditions  found  in  practice.  Instead,  we  consider  the  condition  in  which  no 
heat  sources  are  applied  and  the  temperature  at  the  end  is  held  spatially  constant,  varying 
only  harmonically  in  time.  Then,  rather  than  attempting  to  apply  heat  flux  proportionality 
constants  to  the  test  rig  configurations,  we  need  only  measure  the  temperature  at  the  grip 
location,  scale  the  constant  temperature  solution  appropriately  and  add  it  to  the  solution 
corresponding  to  the  proper  mode  of  heating.  Thus,  only  one  solution  to  account  for 
boundaiy  conduction  needs  to  be  obtained  for  each  frequency. 

In  particular,  we  will  first  investigate  in  detail  the  case  of  symmetrically  heating  a 
specimen  by  heat  lamp  and  cooling  it  by  forced  air.  In  later  chapters,  other  thermal  loading 
cases  will  be  considered  for  comparison.  This  first  method  of  heating  can  be  modeled  by  a 
sinusoidally  varying  uniform  heat  source,  H(x,y,t)  =  ho  cos  tot ,  where  ho  is  spatially 
constant.  So,  in  addition  to  the  constraints  imposed  on  the  stress  and  displacement  vectors 


in  Equations  (1)  through  (1 1),  the  following  conditions  will  be  satisfied  by  the  temperature 
field: 


a.  h  ( x,y  )  =  0 ,  ©o  ( t )  =  cos  tot 

b.  h  ( x.y )  =  ho ,  ©o  ( t )  =  0 


(17) 

(18) 


where  ho  and  Tq  are  constants,  h(x,y)  is  the  heat  source  distribution  of  Equation  (16)  and 
©0  (t)  is  the  edge  temperature  of  Equation  (5). 

With  this  we  have  developed  an  approach  for  establishing  the  effect  of  thermal 
transients  on  the  mechanical  state  of  an  elastic  region  containing  a  crack.  The  temperature 
solutions  are  spatially  one-dimensional  and  are  obtained  in  closed  form.  The  stress  and 
displacement  fields  are  numerically  calculated  with  finite  element  methods.  A  path 
independent  integral  developed  for  thermoelastic  problems  is  used  to  compute  the  stress 
intensity  factor. 


Solutions  to  the  Temperature  Fields 

The  closed  form  representation  of  temperature  for  the  case  with  temperature  varying 
as  TqCos  (ojt)  on  the  boundary  y=L  is  derived  by  standard  separation  of  variables  method 
and  by  invoking  the  Duhamel  integral  [22]: 


^  =  cos  a)t  +  4  ©  ^ 


(-1) 


(n-W 


COS 


nrcy 


Tin 


n= 1,3,5 


(K  + 


sin  ©t  -  ©  cos  ©t) 


(19) 


A  similar  expression  for  the  temperature  in  which  a  uniform  heat  source  varying  as  cos  (©t) 
is  applied  internally  and  the  temperature  at  the  end  is  zero  is  obtained  by  variation  of 
parameters  method  [22]: 
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Figure  2  graphically  depicts  the  temperature  field  according  to  (20)  for  various  fiequencies. 

Solution  to  the  Thermoelastic  Stress  Fields 

The  stress  and  displacement  fields  are  calculated  using  the  finite  clement  method 
which  is  an  indirect  approach  to  solving  partial  differential  equations  describing  a  desired 
quantity  in  a  continuum.  In  this  case,  the  desired  quantity  is  the  stress/displacement  field 
and  the  continuum  is  the  linear  elastic  solid  of  the  two-dimensional  geometry  of  Figure  1. 
In  theory,  the  continuum  is  divided  into  a  number  of  'finite  elements'.  An  approximation 
function  (algebraic  polynomial)  is  chosen  to  represent  the  variation  of  the  desired  quantity 
at  discrete  points  (nodes)  within  the  element  and  on  its  boundary.  The  number  of  nodes  for 
each  element  determines  the  order  of  the  polynomial.  By  using  the  physical  properties  of 
the  continuum  and  the  appropriate  physical  laws,  a  set  of  simultaneous  equations  is 
obtained.  Solution  of  these  equations  results  in  a  representation  of  the  desired  quantity  at 
nodal  locations  within  the  solid. 

The  ABAQUS  finite  element  code  is  a  packaged  program  designed  for  the 
numerical  modeling  of  structural  response.  Its  features  include  its  ability  to  perform  linear 
and  nonlinear  stress  analyses;  stadc  and  dynamic  analyses;  and  coupled  and  uncoupled  heat 
transfer  analyses.  Full  control  of  the  finite  element  selection  is  afforded  the  user.  The  code 
is  capable  of  handling  elements  with  up  to  27  nodes  for  the  three-dimensional  stress 
analysis  case. 

The  uncoupled  thermoelasticity  problem  can  also  be  handled  by  ABAQUS,  making 
it  an  ideal  code  for  the  present  work.  If  the  temperatures  are  specified  at  the  nodal  points 
within  the  material  and  appropriate  material  constants  (including  thermal  expansion 
coefficient)  are  input,  the  code  will  interpolate  the  nodal  temperature  field  to  interior  Gauss 
integration  points  to  calculate  the  thermally  induced  strain  field.  The  stress  field  is  then 
determined  by  the  stress-strain  relations. 

Four-node  bilinear  elements  were  chosen  for  the  analyses.  In  order  to  check  the 
spatial  convergence  of  the  solution,  three  grid  sizes  of  uniform  refinement  were  employed; 
Ax  =  1/6  (coarse  mesh),  1/12  (medium  mesh),  and  1/24  (fine  mesh).  (A  full  description  of 
the  grids  used  in  the  finite  element  analysis  is  contained  in  Appendix  B.)  Square  elements. 
Ax  =  Ay,  were  used  in  all  analyses.  Uniform  refinement  of  the  grid  size  permits  a 
quantitative  evaluation  of  how  well  the  solution  converges.  Moreover,  an  extrapolated 
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value  which  approximates  the  result  of  a  grid  uniformly  finer  than  the  finest  used  in  the 
calculation  can  also  be  computed.  The  extrapolation  scheme  used  models  the  error 
distribution  in  the  numerical  analysis  as 


e 


=  K  -  =  Eo  A*- 


(22) 


where  e  is  the  absolute  percentage  error  in  the  stress  intensity  solution  K  (because  of  the 
singularity  at  the  crack  tip,  this  method  is  inapplicable  for  determination  of  stress  field 
convergence,  rather  it  is  extended  to  the  stress  intensity  solution),  Kext  is  Ae  extrapolated 
solution,  £o  is  the  error  for  the  unit  grid  size,  A  is  the  finite  element  mesh  size,  and  c  is  a 
convergence  measure.  This  extrapolation  scheme  reduces  to  a  Richardson  approximation  if 
c  =1.  Uniform  refinement  implies  the  grid  sizes  are  scaled  proportionally.  In  this  case, 
they  are  reduced  by  a  factor  of  2  with  each  refinement.  Thus  the  extrapolated  value  can  be 
derived  by  algebraic  manipulation  of  the  finite  element  results  from  the  three  grids: 


_  K^Kp.K- 

«‘"K^  +  K,-2K, 


(23) 


Here  the  subscripts  C.  M,  and  F  refer  to  the  coarse,  medium  and  fine  finite  clement 
solutions,  respectively.  The  solution  converges  rapidly  if  c  >1.  Ifc  <0,  the  solution 
diverges.  Between  these  limits,  0  <  c  <1,  impaired  convergence  is  obtained. 

c  =-0.693  ln||*-^| 

(24) 


Obtaining  the  Stress  Intensity  Factor 

The  finite  element  analysis  results  in  solutions  to  the  stress  and  displacement 
vectors  at  discrete  points  within  the  solid.  In  order  to  assess  these  results  in  terms  of  crack 
tip  stress  fields,  stress  intensity  values  are  needed.  The  path  independent  integral  is  a 
method  which  calculates  stress  intensity  values  from  the  data  acquired  in  the  finite  element 
and  finite  difference  analyses.  It  offers  the  distinct  advantage  that  the  errors  in  the  analysis 
stem  only  from  the  numerical  approximations  used  and  not  from  additional  sources  (i.e.  the 
singularity  at  the  crack  tip).  For  this  reason,  path  independent  integrals  offer  a  relatively 
accurate  result  with  only  modest  computational  effoa[13]. 


The  appropriate  integral  for  the  thermoelastic  case  is  developed  following  Sinclair, 
Okajima  and  Griffin  [13],  taking  into  account  the  Duhamel-Neumann  analogy  for  modeling 
thermoelastic  problems  [20].  In  general,  the  analogy  states  that  the  displacements  produced 
by  a  temperature  T  are  the  same  as  those  produced  by  a  body  force  equal  to  -^VT  and 
normal  surface  tractions  equal  to  pT,  acting  on  a  body  of  the  same  shape  but  with  no 
variation  in  temperature.  In  essence,  the  effect  of  the  temperature  change  is  the  same  as  that 
obtained  by  normal  surface  tractions  equal  to  PT  and  body  forces  equal  to  -pVT  .  By 
formulating  the  problem  in  this  way,  it  is  clear  that  we  maintain  the  exact  same 
eigenfunctions  and  eigenvectors  as  those  for  bodies  with  no  temperature  field  (those  of 
Williams  [4])  and  we  need  only  develop  a  new  path  independent  integral  which  takes  into 
account  body  forces  proportional  to  the  gradient  of  the  temperaiure.  To  do  this,  we 
introduce  new  stress  and  displacement  fields,  denoted  by  the  superscript  R  (for  resolved 
tractions  and  body  forces)  and  related  to  those  fields  satisfying  the  equations  of 
thermoelasticity  (7)  through  (1 1)  by  the  relations 


o,  =  -  PT 

Oy  =  Oy*^  -  PT 

T  -  tR 

(25) 


where  P  is  as  defined  in  (8).  The  new  stress  field  satisfies  the  equations  of  linear  elasticity 
in  which  the  body  force  and  normal  surface  tractions  are  supplemented  by  -pVT  and  PT, 
respectively.  These  are  the  same  as  those  presented  earlier  with  the  following  exceptions: 
the  equations  of  equilibrium  in  the  presence  of  a  body  force  field: 

»  *  ■*■  ^xy  >  y  ■  P  T»x  ~  ®  (2g) 

y  +  tfy  ,  ,  -  P  T,y  =0 


on  R,  the  stress-displacement  relations  for  a  homogeneous  and  isotropic,  linear  elastic 
solid. 
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R  ^  M- 

~  (1— 2v)  ^  ^»x  1 

”^xy  ■*■  '^»x) 


(27) 


on  R,  traction  boundary  conditions  on  the  crack  face,  and  on  the  upper  and  right  edges  of 
the  region. 


=  -PT  ,  T^y  =  0  on  y  =  0  (0<x<fl) 
Oy^  =  PT  ,  Tjjy  =  0  on  y  =  L  ( 0  <  x  <  W ) 

ajj^  =  PT  ,  tjjy  =  0  on  x  =  W  (0<y<L) 


(28) 


The  superscripts  on  u  and  v  have  been  omitted  since  the  displacements  compatible  with 
arc  equal  to  those  produced  by  the  original  field  by  the  Duhamel-Neumann  theorem.  We 
now  have  a  problem  we  can  handle  analytically  and,  as  long  as  the  temperature  field  is  not 
more  singular  than  logarithmic  behavior,  there  is  no  contribution  to  the  stress  intensity 
factor  from  the  temperature  field  in  the  superposition  (substituting  Oy  of  Equation  (25)  into 
Equation  (12)).  Thus,  the  stress  intensity  factor  associated  with  the  original  problem  is 
equal  to  that  associated  with  the  fields  and  u.  We  turn  now  to  developing  the  path 
independent  integral  for  regions  containing  body  forces. 

Following  the  procedure  in  Sinclair,  et  al.,[13]  we  invoke  Betti's  reciprocal 
theorem  in  the  plane  for  regions  with  body  forces.  For  the  case  with  a  body  force  field 
equal  to  -pVT,  stress  and  displacement  fields  and  u,  respectively,  Betti's  reciprocal 
theorem  can  be  stated 


u*  + V*  -  a*  u  - v)  n,  +  (Oy*^  V*  +  u*  -  a*  V  -  xjy  u)  Uy  )  dS 


-P  J(T„u*  +T,yV*)dA  =  0 

A 


(29) 


The  starred  functions  are  the  complementary  eigenfunction  stress  and  displacement  fields 
satisfying  the  same  field  equations,  namely  the  complementary  fields  contained  in  Sinclair 
et  al.,  specializing  X=l/2  and  a=7i.  These  fields  are  given  for  this  special  case  in  Appendix 
C.  The  integration  is  performed  in  a  counter-clockwise  direction  along  any  closed  path  F  in 
R.  dS  and  dA  refer  to  line  and  area  elements,  respectively,  nx  and  ny  are  the  components 
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of  the  unit  outward  normal  to  T  and  A  is  the  area  enclosed  by  F.  On  choosing  the  contour, 
we  proceed  as  in  [13]  and  pick  any  path  which  includes  the  crack  faces  and  an  inner 
circular  arc  starting  at  the  crack  face  above  the  axis  of  symmetry  and  ending  on  the  face 
below  it  (Figure  3).  The  choice  of  the  crack  faces  is  based  on  the  homogeneous  traction 
conditions  of  the  complementary  eigenfunction  stress  field  along  the  faces.  We  take  full 
advantage  of  the  crack  face  conditions  by  substituting  the  thermoelastic  stress  field  for 
via  the  relations  of  (25)  since  does  not  generally  equal  zero  on  this  boundary  while  Oy 
is  identically  zero  there.  In  the  limit  as  the  radius  of  the  inner  arc  tends  to  zero,  the  only 
contribution  of  the  inner  arc  to  the  integral  emanates  from  the  singular  parts  of  the  stress 
fields  .  Since  all  stress  fields  now  obey  homogeneous  crack  face  traction  conditions  the 
counter-clockwise  integration  of  the  inner  arc  now  equals  the  integration  along  the  outer 
path  summed  with  the  area  integration.  By  scaling  the  contribution  of  the  circular  arc  so 
that  in  the  limit,  the  stress  intensity  factor  is  recovered,  we  obtain  the  appropriate 
expression  for  K; 


K=|{(a,u*  +  TxyV*  -  oj  u  -  v)  n,  +  (Oy V*  +  x^yU*  -  a*  v  -  x;^  u)  ny  )dS 
-p  J(  T,^  u*  +  T,y  V* )  dA  -i-p  f[(  Tu*)  n^-K Tv*)  Uy]  dS 


(30) 

where  Li  and  Z2  are  the  outer  path  and  crack  face  portions,respectively,  of  the  contour  F . 
This  expression  is  simplified  by  invoking  the  divergence  theorem  to  combine  the  area 
integral  and  the  second  line  integral  as  a  single  area  integral.  Thus,  the  path  independent 
integral  which  defines  the  stress  intensity  factor  for  a  body  subjected  to  a  temperature  T, 
satisfying  the  equations  of  thermoelasticity  (7)  through  (11),  is 


K=J{(a,u*  +  x^v*  -  -  x;^  v)  n,+  (OyV*  +  x^u*  -  a*  v  -  x;jy  u)  ny  )dS 

-hPJt(u*.*  +  v*,y)dA 


(31) 


The  resulting  path  independent  integral  is  the  sum  of  an  area  integral  containing  the 
temperature  field  and  a  line  integral  containing  the  stress  and  displacement  fields. 
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Numerical  computation  of  K  requires  compatible  quadratures  to  accurately  determine  the 
stress  intensity  factor.  Both  the  temperature  and  displacement  fields  are  determined  at 
nodal  locations,  while  the  stress  field  is  most  accurate  at  the  interior  Gauss  integration 
points.  However,  the  ABAQUS  code  does  interpolate  these  stress  values  to  node  locations 
if  so  desired  and  integration  was  performed  using  these  interpolated  nodal  values. 

Integration  using  Simpson's  rule  is  based  on  the  use  of  parabolic  arcs  rather  than 
straight  lines  to  approximate  the  integrand  [23]  and  is  accurate  to  order  (Ax)^  .  This 
technique  uses  nodal  values  of  the  function.  Thus,  it  is  well  suited  for  the  evaluation  of  the 
line  integral  since  this  approach  does  not  require  evaluation  of  the  function  near  its 
singularity,  the  crack  tip. 

Solution  of  the  area  integral  is  not  as  straightforward.  Upon  evaluation  of  uij*  ,  the 
area  integral  becomes  (in  cylindrical  coordinates  with  the  origin  at  the  crack  tip) 


f  P  T  cos  dr  d9 

i  ^ 


(32) 


for  the  Mode  I  stress  intensity  factor.  The  behavior  of  the  complementary  eigenfunction 
displacement  field  results  in  an  integrable  square  root  singularity  in  the  integrand  in  one 
coordinate  (the  other  remains  nonsingular).  Because  of  this  singularity  and  the  fact  that  the 
integrand  must  be  numerically  evaluated  throughout  the  area,  a  quadrature  which  can 
handle  this  singularity  is  required.  A  Gauss-Legendre  quadrature  was  developed  following 
[24]  to  approximate  the  general  form  of  the  integral: 


(33) 


Defining  the  integral  in  this  manner  removes  the  difficulty  in  integrating  the  singularity 
since  we  require  f(x)  to  be  regular.  If  f(x)  is  a  polynomial,  this  approximation  is  exact 
when  used  with  the  proper  weights,  wj ,  and  abscissas,  Xj .  For  this  particular  square  root 
weight  function,  the  abscissas  are  the  zeros  of  the  Legendre  polynomials  ,  Pn  (^)  [25]. 
Values  for  the  weights  and  abscissas  were  calculated  for  up  to  96  interior  points  and  are 
tabulated  in  Appendix  D. 

For  the  nonsingular  integration  (over  6) ,  standard  Gauss-Legendre  quadrature  is 
used.  Tabulated  values  of  weights  and  abscissas  for  this  method  are  common.  Numerical 
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results  of  the  area  integral  using  these  two  Gauss-type  quadratures  converged  rapidly  when 
tested  on  several  problems: 

f(x)=  1.0 
f(x)  =  1  -  2x 
f(x)  =  5  -  6x^ 
f(x)  =  Tx 
f(x)  =  In  (x) 

With  the  first  three  functions,  the  value  converged  to  the  sixth  decimal  point  with  only  six 
interior  Gauss  points.  Convergence  on  the  last  two  functions  was  somewhat  slower.  For 
the  square  root  function,  the  value  converged  to  the  fourth  decimal  point  after  48  interior 
points;  for  the  logarithmic  function,  the  value  converged  to  the  third  poin  after  48  interior 
points. 

In  order  to  verify  the  stress  analysis  and  stress  intensity  calculation  techniques,  we 
perform  several  test  problems.  In  each  of  the  test  problems,  static  loadings  are  applied  to 
center-cracked  geometries.  The  loading  conditions  for  each  of  these  problems  are  detailed 
in  Table  1  and  include  (i)  uniaxial  tension  applied  at  y=L,  no  temperature  distribution,  (ii) 
displacements  applied  on  the  boundary  y«L,  no  temperature  distribution,  and  (iii)  a 
spatially  constant  temperature  distribution  with  the  ends  fixed  from  vertical  displacement. 
Normalized  stress  intensity  factors,  taken  from  the  literature  are  also  presented  in  Table  1 
and  are  exact  for  infinitely  long  specimens.  Note  that  width  correction  factors,  taken  from 
the  literature,  are  also  included  in  the  table. 

In  general,  excellent  agreement  with  the  quoted  values  for  all  test  problems  is 
observed  in  the  results  which  are  listed  in  Table  2.  The  calculations  were  performed  in 
plane-strain  with  a  Poisson  ratio  of  0.3.  The  results  are  tabulated  as  normalized  stress 
intensity  factors  divided  by  the  appropriate  reported  solution  described  in  Table  1  such  that 
an  exact  value  equals  unity.  The  normalized  crack  length,  a/W  =  1/3  for  all  cases  listed  in 
Table  2.  The  normalizing  factor  is 

Kj 

Co-Tm 


K  = 


(34) 
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where  Go  refers  to  the  applied  load.  For  test  case  (i),  Oq  is  simply  the  traction  applied  at 
y=L.  For  test  case  (ii),  Oq  is  given  by 


1-v^  ^  (35) 

For  test  case  (iii),  Gq  is  given  by 
gET^ 

1-v  (36) 

The  normalized  crack  length,  a/W  ,  for  the  values  quoted  is  1/3.  In  addition,  stress 
intensity  values  were  obtained  for  o/W  =  1/2  and  the  accuracy  of  the  results  was  found  to 
be  consistent  with  that  computed  for  the  shorter  crack  length. 

In  addition  to  the  individual  values  generated  by  coarse,  medium  and  fme  grids,  the 
extrapolated  stress  intensity  factors  obtained  by  using  Equation  (23),  and  convergence 
measures  from  (24)  are  included  in  the  tables.  Further,  path  independence  of  the  integral 
was  demonstrated  by  performing  the  integration  along  two  different  paths:  one  which  takes 
advantage  of  one  boundary  of  the  specimen,  another  which  remains  in  specimen  interior. 
The  paths  of  integration  are  illustrated  in  Figure  4.  In  each  case,  the  two  contours  converge 
to  similar  values  of  the  stress  intensity  factor.  The  variation  of  the  two  paths  is  of  the  order 
of  20  per  cent  for  the  coarse  grid  and  converges  to  approximately  2  per  cent  for  the  fine 
grid  (less  than  1  percent  for  the  extrapolated  values).  However,  slower  convergence  was 
observed  with  the  contour  on  the  boundary,  probably  due  to  the  difficulty  in  satisfying  the 
boundary  conditions  in  the  finite  element  analysis,  resulting  in  relatively  large  errors  along 
the  portion  of  the  integration  path  which  is  on  the  boundary.  For  this  reason,  the  interior 
contour  was  selected  for  all  successive  calculations. 

Numerical  Results 

Keeping  the  normalized  crack  length  constant  and  equal  to  1/3,  dimensionless  stress 
intensity  factors  were  computed  for  dimensionless  frequencies  ranging  from  0  to  8.2  This 
range  is  representative  of  typical  thermal  cycle  frequencies  used  in  fatigue  testing.  For 
example,  the  high-strength  nickel  alloy  IN  100  has  a  diffusivity  of  0.00649  in.2/s3.  If  the 


^The  dimensionless  frequency,  to*,  is  related  to  the  frequency,  o),  by  the  relation  to*  « to  /  (D  n) . 
3For  k»100  Btu  in/(ft2  hr  ©F).  p.0.28  Ib/in^.,  Cp=0.106  cal/g®C  and  D=k/(p  Cp). 
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total  specimen  length  is  2  inches  (2L  =  2  in.),  the  normalized  frequency  a)*=2.5 
corresponds  to  a  real  frequency  of  approximately  o)=0.5  cycles/min.  Results  of  the 
computations  are  tabulated  in  Table  3  for  the  case  with  the  sinusoidally  varying  end 
temperature  and  in  Table  4  for  the  case  with  sinusoidally  varying  uniform  heat  source.  As 
listed,  the  results  are  the  dimensionless  magnitude,  K,  and  phase  lag,  d>,  of  the  harmonic 
stress  intensity  factor; 


K, 


(i-v) 


Vtw 


=  K  cos  (cot  -  O) 


(37) 


where  the  dimensioning  parameter  Xq  refers  to  the  end  temperature  Tq  in  Table  3  and  the 
magnitude  of  the  heat  source  hoL^  in  Table  4.  If  Kj  and  Kc  are  the  dimensionless  stress 
intensity  factors  associated  with  the  sin  (cot)  and  cos  (tot)  components  of  the  temperature 
field  (Equations  (14)  and  (15)),  then  the  dimensionless  magnitude  K  is  obtained  by  taking 
the  square  root  of  the  sum  of  the  squares  of  Ks  and  IQ: 


k  =  VkJ+iq^ 


(38) 


and  the  phase  lag,  O,  is  obtained  thus; 


A  K, 

<P  =  arctan-^ 
Kc 


(39) 


Because  the  temperamre  in  the  vicinity  of  the  crack  tip  is  of  interest.  Tables  3  and  4 
also  include  the  magnitude,  Ti,  and  phase,  <l>i,  of  the  dimensionless  crack  tip  temperature: 

—  =  Tj  cos  (tot  -  Oj) 


(40) 
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In  all  cases,  the  values  quoted  are  obtained  from  the  numerical  analyses  using  the 
extrapolation  method  described  previously  (23).  The  convergence  of  the  stress  intensity 
factors  was  found  to  be  well  behaved. 

It  is  of  interest  to  assess  the  variation  in  stress  intensity  with  frequency  for  the 
uniform  heating  condition  in  which  the  stress  intensity  factor  is  normalized  by  the 
temperature  at  the  centerline  (the  crack  tip).  With  this  normalization  we  can  observe  the 
trend  in  stress  intensity  with  frequency  when  a  specific  temperature  is  maintained  within  the 
sample  if  the  temperature  at  the  end  is  maintained  as  zero.  These  data  are  plotted  in  Figure 
5.  It  is  noted  that  the  maximum  stress  intensity  occurs  for  the  static  case,  O)  =  0 . 

Using  linear  superposition  we  can  evaluate  the  Mode  I  stress  intensity  for  different 
fatigue  conditions  in  which  uniform,  symmetric,  heating  is  applied.  By  measuring  the 
temperature  and  phase  lag  at  two  locations  within  the  sample  (for  example,  at  the  crack  tip 
and  at  a  distance  L  from  it),  we  can  derive  the  stress  intensity  factor  for  specific  frequencies 
of  thermal  cycle  using  the  data  presented  in  Tables  3  and  4.  First  consider  a  specimen  that 
is  uniformly  heated  and  in  which  the  temperature  at  the  crack  tip,  Ti  is  measured  with 
respect  to  the  end  temperature,  Tq.  The  phase  lag  of  the  measured  crack  tip  temperature  is 
then  called  d)i.  With  this  information,  we  can  combine  the  heat  source  results  with  the 
boundary  temperature  results  to  estimate  the  contribution  of  thermal  cycling  on  stress 
intensity.  To  maintain  the  crack  tip  temperature  at  Ti,  a  sufficient  magnitude  of  the  heat  is 
required  at  a  specific  phase  lag  from  the  boundaiy  temperature.  Since  we  have  analyzed  the 
temperatures  at  the  crack  tip  in  all  the  problems,  we  need  only  find  the  magnitude  of  the 
heat  source  and  the  phase  required  to  maintain  the  crack  tip  at  a  specific  temperature.  In 
other  words,  if  Tib  and  Oib  are  the  magnitude  and  phase  of  the  crack  tip  temperature  in 
the  case  where  T©  is  specified  on  the  boundary  y=L,  and  Tih  and  Oih  are  the  analogous 
values  in  the  cases  where  a  symmetric  heat  source  is  present,  then 

Tj  cos  (cot  -  <Dj )  =  Tjg  cos  ( cot  -  0,g )  +  y Tj„  cos  ( cot  -  -  k )  ... 


where  y  is  a  scaling  factor  for  the  heat  source  center  temperature  and  k  is  the  required 
phase  lag  of  the  heat  source.  These  two  parameters  are  determined  by  the  measurements  Ti 
and  <I>i,  and  the  results  presented  in  Tables  3  and  4: 


(42) 
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and 


K  =  arctan 


Ti  sin  Oj-TjpSin 
Tj  cos  Oj-Tjg  cosOjg 


(43) 


By  establishing  Y  and  K,  the  stress  intensities  may  also  be  found.  If  Ktot  represents 
the  total  dimensionless  stress  intensity  factor  in  the  specimen  with  crack  tip  temperature  Ti 
cos  (cot  -  Oi)  relative  to  the  end  temperature,  T©  cos  cot,  then 


=  V  (YKh)'  +  2YK3  K„  cos  (®3  -  C>„  -K ) 


(44) 


where  Kb  and  <I>b  correspond  to  the  magnitude  and  phase  of  the  stress  intensity  with 
boundary  temperature  prescribed  and  Kh  and  Oh  are  the  respective  values  of  the  stress 
intensity  with  an  applied  heat  source.  For  the  static  case,  if  the  crack  tip  temperature  is 
maintained  at  the  same  level  as  the  end  temperature  then  the  stress  intensity  is  zero  since  the 
temperature  is  uniform  within  the  plate.  This  is  consistent  with  Equations  (41)  through 
(44).  If  Ti  =  To  and  00  =  0  (Oi  =0),  then  from  (42),  y  =  0.  Since  Kb  is  0  in  the  static 
case,  we  have  Ktot  =  (J* 

Two  examples  of  this  type  of  superposition  are  considered  and  results  are  plotted  in 
Figures  6  and  7.  In  these  plots,  the  combined  stress  intensity  ,  Ktot.  is  plotted  as  a 
function  of  T i/To  for  several  frequencies.  Two  phase  lags  are  illustrated  and  correspond  to 
the  crack  tip  temperature  in  phase  and  ISOo  out  of  phase  with  the  end  temperature.  It  is 
noted  that  the  stress  intensity  is  zero  for  all  frequencies  in  the  uniform  heating  case  when 
the  crack  tip  temperature  is  equal  to  and  in  phase  with  the  end  temperature.  What  this  says 
is  that  all  temperatures  in  the  plate  vary  according  to  Tq  cos  cot  when  this  constraint  is 
imposed  at  both  the  end  of  the  specimen  and  at  the  crack  tip.  In  order  for  this  to  occur 
numerically,  the  sum  of  the  phases  Oh.  tmd  K  must  add  to  ±  180°  in  Equation  (44). 
(This  is  true  as  long  as  y  is  non-zero  which  is  the  case  for  nonzero  frequency. )  The  sum 
of  Oh  and  Ob  from  Tables  3  and  4  is  equal  to  -90°  to  the  third  decimal  point  for  each 
nonzero  frequency.  Determining  k  from  Equation  (43)  is  straightforward  and  equal  to  -90° 
leading  to  a  sum  of  -180°.  These  data  confirm  the  accuracy  of  the  numerical  values  in  the 
tables. 


Summary  of  Results 

In  this  section,  we  have  developed  a  method  to  estimate  the  variation  of  thermal 
fatigue  stress  intensities  as  a  function  of  frequency  for  the  basic  problem  of  a  uniform  heat 
source  in  a  center-cracked  plate.  Resulting  values  of  stress  intensity  were  computed  for  a 
range  of  frequencies  and  show  that  the  static  case,  O)  =  0,  gives  the  strongest  stress 
intensity.  Further,  the  stress  intensity  depends  on  the  phase  and  magnitude  of  the 
temperature  at  the  specimen  end  relative  to  the  temperature  along  the  centerline.  In  order  to 
estimate  the  stress  intensity  for  specific  fatigue  test  cases,  laboratory  data,  specifically  the 
relative  temperature  values,  are  required. 
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Table  1 

Loading  conditions  and  reponed  stress  intensity  solutions  for  test  problems 
The  plane  strain  solutions  are  normalized  such  that 


Test 

Case 

Applied 

Conditions 

Reported  Solutions 
(Center  crack  geometry) 

Reference 

(i) 

Cy  =aoaty=L 

K=Fi(a/W) 

H 

II 

o 

Fi(l/3)=1.07 

[26] 

(ii) 

v=Vo  at  y=L 

T  =  0 

Vo  E 

q 

II 

[26] 

(iii) 

v=0  at  y=L 

T  =  To 
oETo 

^0-  i_v 

II 

o 

[26] 
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Table  2 

Results  of  test  problem  calculations  using  center-crack  geometry 
Reported  as  normalized  stress  intensity  factor,  precise  value  for  Kext  is  1.0 
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Table  3 

Dimensionless  stress  intensity  f£u:tor  and  temperature  field  results  for 

center-crack  geometry, 
end  temperature  prescribed  non-rero,Xo  =  Tq 


to* 

K 

Ti 

Ol 

(X  10-3) 

(radians) 

(radians) 

0 

0.0 

-1.571 

1.0 

0.0 

0.5 

47.623 

-0.937 

0.841 

0.716 

1.0 

69.175 

-0.528 

0.610 

1.202 

1.5 

76.899 

-0.277 

0.452 

1.532 

2.0 

79.348 

-0.102 

0.349 

1.784 

3.0 

78.512 

0.139 

0.229 

2.183 

4.0 

74.793 

0.312 

0.163 

2.513 

8.0 

55.627 

0.710 

0.0577 

3.544 
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Table  4 

Dimensionless  stress  intensity  factor  and  temperature  field  results 
for  center  crack  geometry,  uniform  heat  source 

to  =  hoL2 


CO* 

K 

O 

Ti 

‘I>1 

(X  10-3) 

(radians) 

(radians) 

0 

36.125 

-3.141 

0.5 

0.0 

0.5 

30.391 

-2.505 

0.421 

0.584 

1.0 

22.071 

-2.098 

0.307 

0.940 

1.5 

16.359 

-1.847 

0.229 

1.140 

2.0 

12.657 

-1.673 

0.179 

1.447 

3.0 

8.358 

-1.432 

0.122 

1.407 

4.0 

5.980 

-1.260 

0.0904 

1.486 

8.0 

2.416 

-0.864 

0.0419 

1.592 
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Figure  2.  Magnitude  of  the  temperature  field  induced  by  uniform  heating^. 


^The  dimensionless  frequency,  co*.  is  related  to  the  fiequency.  u,  by  the  relation  o)*  *  (oL^/  (Dn). 


a.  Interior  Path 


b.  Boundary  Path 


Figure  4.  Integration  paths  for  the  test  problems.  Areas  shown  represent  the 
upper  quadrant  for  the  center  crack  geometry. 


32 


0  2  4  6  8 

Normalized  frequency 

Figure  5.  Amplitude  of  the  stress  intensity  normalized  by  the  crack  tip  temperature 

for  the  uniform  heating  case. 


Figure  6.  Amplitude  of  the  stress  intensity  as  a  function  of  the  crack  tip  temperature.  The 
temperature  at  the  crack  tip,  Ti,  is  maintained  by  uniform  heating  and  is  in  phase 
with  the  end  temperature,  Tq. 
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T,/T, 


Figure  7.  Amplitude  of  the  stress  intensity  as  a  function  of  the  crack  tip  temperature  and 
dimensionless  frequency.  The  temperature  at  the  crack  dp,  Ti,  is  maintained  by 
uniform  heating  and  is  180°  out  of  phase  with  the  end  temperature,  Tq. 


ffl.  EXTENSIONS  OF  THE  UNIFORM  HEAT  SOURCE  PROBLEM 


In  this  section,  we  consider  several  extensions  to  the  uniform  heat  source  problem 
formulated  and  solved  in  Chapter  n.  First,  because  fatigue  test  specimens  are  not  all 
geometrically  equivalent  to  the  center-crack  geometry,  we  will  investigate  some  geometrical 
differences:  edge-crack  geometries  and  plate  and  crack  length  scaling  dependence,  in 
particular.  Secondly,  because  specimen  heating  may  not  always  be  symmetric,  we  will 
consider  antisymmetric  temperature  fields  which  give  rise  to  Mode  n  stress  intensities. 
These,  when  combined  linearly  with  the  symmetric  heating  results,  will  lead  to  predictions 
of  the  stress  intensity  in  more  general  thermal  fatigue  situations.  Differences  in  problem 
formulation  and  solution  will  be  discussed  in  the  pertinent  sections,  below. 

Geometrical  variations 
a.  Edge-Crack  Geometry 

In  fatigue  testing,  cracks  can  be  located  on  the  specimen  edge  or  in  its  center.  The 
important  distinction  in  the  crack  location  is  that  the  edge-cracked  component  corresponds 
to  a  simply-connected  region  while  the  center-cracked  component  corresponds  to  a 
multiply-connected  region.  In  our  investigation  of  this  geometric  variation,  we  use  the 
same  formulation  as  that  already  presented  and  the  edge-crack  geometry  of  Figure  8.  The 
region  R  of  the  analysis  is,  by  symmetry,  the  upper  half  of  the  rectangular  plate  shaded  in 
the  figure  and  corresponds  to  the  same  region  investigated  in  the  center-crack  problem: 

R  ={x,y  l0SxSW,0^y^L} 

A  crack  of  length  a  is  contained  on  the  boundary  0  ^  x  ^  a  along  the  line  y  =  0. 
The  boundary  conditions  along  the  edge  containing  the  crack  are  the  only  variation  from  the 
formulation  already  presented.  The  condition  on  the  temperature  problem  is  that  the  edge 
remains  insulated  which  is  identical  to  the  condition  of  symmetry  imposed  on  the  center- 
crack  problem  (4).  Because  of  this,  the  temperature  fields  for  the  two  geometries  are 
identical  and  will  not  be  reiterated  here. 
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The  conditions  for  the  stress  field,  however,  are  different  along  the  edge  containing 
the  crack.  On  this  boundary,  traction-&ee  conditions  are  imposed; 

Ox  =  0,  Xxy  =  0  onx  =  0  (0<y<L) 

The  remaining  conditions.  Equations  (9)  and  (1 1),  still  hold. 

With  this,  we  have  the  formulation  for  the  edge-cracked  geometry.  Solution  is 
accomplished  with  the  same  methods  described  in  the  previous  section:  temperature  values 
at  nodal  locations  are  calculated  from  the  closed  form  temperature  fields  of  Equations  (19) 
and  (20)  and  are  used  as  input  to  the  finite  element  calculations.  The  path  independent 
integral  derived  in  Equation  (31)  is  again  used  to  calculate  the  dimensionless  stress  intensity 
factors. 

To  verify  the  stress  intensity  calculations  for  the  edge-cracked  geometry  several  test 
problems  are  again  performed.  The  three  static  test  problems  described  for  the  center-crack 
geometry  are  again  used  in  addition  to  a  fourth  test  case  in  which  is  applied  a  static 
temperature  distribution  proportional  to  x/W  with  ends  fixed  from  vertical  displacement. 
Table  5  contains  descriptions  of  the  loading  conditions  for  each  of  the  test  problems  and 
width  and  bending  correction  factors  taken  from  the  literature.  Results  of  the  test  cases, 
normalized  such  that  an  exact  value  equals  unity,  are  listed  in  Table  6  for  coarse,  medium 
and  fine  mesh  sizes.  The  calculations  were  performed  in  plane-strain  with  Poisson's  ratio 
equal  to  0.3.  Unless  specified,  results  quoted  are  for  normalized  crack  lengths  of  1/3.  In 
general,  excellent  agreement  with  the  quoted  values  for  aU  test  problems  is  observed  for  the 
frne  mesh,  with  errors  of  less  than  two  per  cent  generally  observed. 

Dimensionless  stress  intensity  factors  for  the  posed  edge-crack  problem  were 
computed  for  the  same  range  of  frequencies  as  the  center-crack  problem  and  are  tabulated  in 
Table  7  for  the  case  with  sinusoidally  varying  end  temperature  and  in  Table  8  for  the  case 
with  sinusoidally  varying  uniform  heat  source.  Again,  the  listed  values  are  the 
dimensionless  magnitude,  K,  and  the  phase  lag,  O,  of  the  harmonic  stress  intensity  factor. 
It  is  noted  that  the  stress  intensity  values  obtained  are  much  lower  than  those  of  the  center- 
cracked  geometry,  often  by  a  factor  of  five.  Similar  trends  of  decreasing  stress  intensity 
beyond  a  dimensionless  frequency  of  about  2.0  are  observed. 

b.  Oack  and  sample  length  scaling  effects. 
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Although  nondimensionalizing  provides  a  way  to  generalize  the  problem  so  that 
many  different  sizes  of  samples  can  be  evaluated  with  only  one  calculation,  it  does  not  take 
into  account  variations  in  ratios  of  crack  length  and  the  sample  dimensions.  In  this  section, 
we  investigate  the  dependence  of  length  scaling  on  the  stress  intensity  factor  by  computing 
the  stress  intensity  factor  as  a  function  of  frequency  for  shorter  relative  crack  lengths, 
fl/W=l/6,  and  longer  relative  sample  lengths,  L/W  =  4.0.  For  comparison,  the  case  in 
which  the  heat  source  varies  sinusoidally  was  used  in  these  scaling  calculations  .  Results 
obtained  by  setting  the  normalized  crack  length  equal  to  half  that  used  in  the  original 
calculations,  are  tabulated  in  Table  9.  Results  of  doubling  the  relative  sample  length  are 
listed  in  Table  10. 

The  stress  intensity  results  listed  are  nondimensionalized  by  the  square  root  of  the 
crack  length.  From  Tables  4  and  9,  since  the  values  of  normalized  stress  intensity  are 
nearly  the  same,  it  is  apparent  that  the  stress  intensity  factor  scales  approximately  as 
(the  crack  length  normalization  factor).  Further,  this  scaling  is  independent  of  the 
fiequency.  The  scaling  is  consistent  with  stress  intensities  generated  by  mechanical 
loads. 

To  assess  the  effects  of  scaling  the  sample  length,  we  compare  the  results  of  Tables 
4  and  10.  Doubling  the  length  of  the  sample  relative  to  the  width  apparently  reduces  the 
stress  intensity  in  the  sample  by  a  factor  of  approximately  four,  (L/W)2,  for  the  same 
temperature  value  at  the  crack  tip.  It  is  noted,  however,  that  this  solution  may  be  particular 
to  this  situation  and  may  not  be  a  general  result  since,  by  lengthening  the  sample,  the 
temperature  constraint  at  the  end  y=L  would  probably  be  affected  as  well. 

Asymmetrical  temperature  distributions 

The  equations  and  boundary  conditions  which  define  the  temperature  fields  are 
linear  and  so  the  fields  themselves  may  be  combined  via  linear  superposition.  This  implies 
that  an  asymmetric  temperature  field  can  be  decomposed  as  the  sum  of  a  symmetric  field 
and  an  antisymmetric  field.  If  we  consider  the  uniform  heat  source  and  the  symmetric  end 
temperature  conditions  as  the  symmetric  contribution  and  a  linearly  varying  heat  source  and 
antisymmetric  end  temperature  conditions  as  the  antisymmetric  contributions,  we  can  obtain 
the  coefficients  to  the  Hrst  terms  in  a  Taylor  series  approximation  to  a  more  general  class  of 
problems.  These  antisymmetric  fields  will  be  developed  and  investigated  in  this  section. 

The  constraints  imposed  in  the  formulation  of  the  antisymmetric  problem  are  that 
the  temperature  at  the  end  y=+L  is  the  negative  of  the  temperature  at  the  end  y»-L.  This 
gives  the  antisymmetric  condition  of  temperature  equal  to  zero  along  the  x^axis,  not 
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including  the  crack  (a<x<  W).  Along  the  crack  faces,  the  condition  of  insulation  (nonnal 
gradient  equal  to  zero)  still  holds.  Thus,  a  discontinuity  in  the  temperature  boundary 
conditions  occurs  at  the  point  x  =  a,  and  the  temperature  fields  cannot  be  obtained  in  closed 
form.  In  addition,  the  heat  source  will  be  imposed  as  one  that  varies  linearly  in  the  y- 
direction: 


h  (x,y)  =  ho  |- 


(47) 


for  a  constant  ho.  The  coupled  equations  defining  the  temperature  field  (16)  are  elliptic  and 
so  the  fields  Ts  and  Tc  may  be  calculated  using  finite  difference  methods  as  described 
below. 

The  boundary  conditions  pertaining  to  the  stress  field  are  identical  to  those  of  the 
symmetric  problem  with  the  exception  of  those  along  the  x-axis  in  the  region  a  <  x  <  W. 
Again,  antisymmetric  conditions  are  imposed : 

Oy  =0  and  u  =  0ony  =  0  (a<x<W) 

The  stress  field  solution  will  again  be  obtained  numerically  using  the  finite  element  method 
already  discussed. 

a.  The  finite  difference  method 

Finite  difference  methods  employ  Taylor  series  expansions  to  express  derivatives  of 
analytic  functions  by  arithmetic  operations.  In  our  analysis  we  employ  representations  of 
derivatives  involving  error  of  order  (Ax)2  When  these  difference  formulas  are  substituted 
into  the  coupled  equations  for  temperature,  a  set  of  linear  algebraic  equations  is  obtained. 
The  method  originally  chosen  for  solving  these  equations  was  Gauss-Seidel  iteration  with 
an  absolute  convergence  criterion.  The  iterative  method  starts  with  a  guess  for  the 
unknown  vector  and  cycles  through  the  equations  replacing  the  solution  for  the  unknowns 
until  each  solution  satisfies  the  convergence  criterion.  This  continues  until  a  change  in  each 
unknown  from  the  previous  iteration  to  the  current  iteration  is  less  than  an  acceptable  value. 
This  value  was  set  to  be  no  greater  than  10"®.  With  a  criterion  this  small,  however, 
convergence  was  too  slow  for  the  two-dimensional  fields  and  with  greater  values,  errors  in 
the  temperatures  were  unacceptably  large.  Because  of  this,  an  alternate  method,  Gauss- 
Jordan  elimination,  was  selected  and  used  for  the  more  complex  temperature  fields.  This 
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technique  uses  row  and  column  combination  until  a  lower-  or  upper-triangular  array  of 
equations  results.  The  unknown  vector  is  then  solved  for  directly. 

Using  Gauss-Jordan  elimination  provides  a  further  advantage  over  iteration  in  that 
accuracy  of  the  solution  depends  only  on  the  grid  size  used  and  round-off  errors  rather  than 
the  additional  error  associated  with  iteratively  solving  the  set  of  algebraic  equations  with  too 
large  a  convergence  criterion.  In  order  to  check  the  spatial  convergence  of  the  solution, 
three  grid  sizes  of  uniform  refinement  were  employed  in  the  analysis:  Ax  =  1/6  (coarse 
mesh),  1/12  (medium  mesh),  and  1/24  (fine  mesh).  Square  elements.  Ax  =  Ay,  were  used 
in  all  temperature  analyses.  Uniform  refinement  of  the  grid  size  permits  use  of  the 
extrapolation  scheme  discussed  in  the  finite  element  section  to  model  the  error  distribution 
in  the  numerical  analysis  as 


e 


(49) 


where  e  is  the  error  in  the  temperature  solution  T  at  a  particular  node  location.  Text  is  the 
extrapolated  solution,  eo  is  the  error  for  the  unit  grid  size,  A  is  the  finite  difference  mesh 
size,  and  c  is  the  convergence  measure. 

To  verify  the  numerical  solution  method,  we  use  several  techniques.  First,  the 
closed  form  one-dimensional  representations  of  the  symmetric  temperature  problem 
(Equations  (19)  and  (20))  are  compared  with  the  values  obtained  inputting  the  appropriate 
boundary  conditions  and  forcing  function  for  the  synunetric  problems.  Secondly,  a  one¬ 
dimensional  antisymmetric  problem  in  which  the  conditions  T=0  at  y=0  and  T=To  at  y=L 
was  formulated  and  solved  in  closed  form  and  compared  with  a  numerical  solution. 
Finally,  boundary  values  for  three  different  known  temperature  distributions  were  input  to 
the  code  and  the  resulting  numerical  solutions  compared  with  the  known  distributions. 
These  three  temperature  disuibutions  and  their  respective  forcing  functions  are: 

T  (x.y)  _  ,  x^-f  y^  .  h  (x,y)  _  , 

hoU  hoL 

T(x,y)  _  .  y'*  .  h  (x,y)  _  -  i 
t.T2  "  4’  .,2 


(51) 
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T  (r,  6) 
To 


Vr  sin  y  ;  h  (r.O)  =  0 


(52) 


where  the  spatial  variables  x  and  y  are  normalized  to  the  specimen  length,  L,  r  =  (x2+y2)l/2 
and  6=aictan  (x/y).  In  general,  all  three  verification  techniques  yielded  numerical  solutions 
which  converged  rapidly  to  the  desired  solutions,  within  one  per  cent  error  for  the  mesh 
sizes  1/6, 1/12  and  1/24. 


b.  Numerical  results 

Again  the  crack  length  was  maintained  as  W/3  and  dimensionless  Mode  11  stress 
intensity  factors  were  computed  for  dimensionless  frequencies  ranging  from  0  to  8  with  the 
same  dimensioning  parameters  as  those  in  the  symmetric  problem: 


Kn 


-gETp 

(l-v) 


/m 


K  cos  (tot  -  O) 


(53) 


Results  of  the  computations  are  tabulated  in  Table  11  for  the  case  with  sinusoidally  varying 
end  temperature  and  in  Table  12  for  the  case  with  sinusoidally  varying  linear  heat  source^. 
The  tables  also  contain  the  temperature  magnitudes  and  phase  lags  (Ti  and  <I>i)  at  the 
location  x  =  0,  y  =  and  not  at  the  crack  tip  since  the  temperature  is  constrained  to  be 
zero  there. 

In  general,  the  stress  intensity  magnitudes  vary  much  less  with  frequency  than 
those  of  the  symmetric  problem  and  are  generally  much  lower  at  specific  frequencies  than 
the  corresponding  Mode  I  stress  intensities. 


Summary  of  Results 

In  this  section  we  considered  several  extensions  to  the  basic  problem  formulated 
and  analyzed  in  the  previous  section.  First,  in  considering  edge-crack  geometries,  we 
found  that  the  edge-crack  region  results  in  stress  intensities  much  lower  than  the  center- 
crack  region.  Secondly,  we  found  that  the  stress  intensity  is  proportional  to  a^/^.  This  is 
consistent  with  stress  intensities  induced  by  mechanical  loads.  In  a  particular  problem  in 
which  specimen  length  scaling  was  examined,  the  stress  intensity  scaled  as  the  inverse  of 


^The  temperature  fields  for  these  specific  problems  were  found  to  converge  adequately. 
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the  square  of  the  length.  However,  this  is  not  considered  a  general  result  since  the 
temperature  constraint  at  the  specimen  end  would  probably  change  with  the  scaling. 
Finally,  antisymmetric  temperature  fields  were  evaluated  and  Mode  n  stress  intensities 
calculated.  In  general,  we  found  these  stress  intensities  to  be  smaller  than  the  Mode  I  stress 
intensities  but  in  order  to  evaluate  the  significance  of  the  Mode  II  values,  representative 
cases  would  have  to  be  evaluated. 
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Tables 

Loading  conditions  and  reported  stress  intensity  solutions  for  test  problems 
The  plane  strain  solutions  are  normalized  such  that 


Test 

Case 

Applied 

Conditions 

Reported  Solutions 

(Edge  crack  geometry) 

Reference 

(i) 

Oy  —  (Jq  at  y— L 

K=F2(aAV) 

T  =  0 

F2(1/3)=1.7864 

[26] 

(ii) 

v=Vo  at  y=L 

K=G(o/WJL/W) 

T  =  0 

L 

G(l/3,2)=1.2467 

[27] 

(iii) 

v=0  at  y=L 

K=G(flAVL/W) 

T  =  To 
aETo 

^0  -  i_v 

G(1/3,2)=1.2467 

[27] 

(iv) 

v=0  at  y=L 

T  =  2T  — 

°  W 

aETo 

^0--  i_v 

K=0.5 

[15] 
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Table  6 

Results  of  test  problem  calculations  using  edge-crack  geometry 
Reported  as  normalized  stress  intensity  factor,  precise  value  is  1.0 


Test 

Case 

Coarse 

Grid 

Medium 

Grid 

Fine 

Grid 

Kext 

c 

(i) 

0.8000 

0.9403 

0.9660 

0.9718 

1.17 

(ii) 

0.8711 

0.9808 

0.9839 

0.9840 

2.47 

(iii) 

0.8510 

0.9843 

0.9849 

0.9850 

3.64 

(iv)7 

0.9280 

0.9746 

0.9814 

0.9826 

1.33 
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Tabic? 

Dimensionless  stress  intensity  factor  results 
for  the  sinusoidally  varying  end  temperature, 
Edge  crack  geometry 
'Co  =  To 


CD* 

K 

O 

(X  10-3) 

(radians) 

0 

0.0 

-1.571 

0.5 

9.324 

-0.883 

1.0 

13.504 

-0.415 

1.5 

15.060 

-0.107 

2.0 

15.510 

0.125 

3.0 

15.252 

0.481 

4.0 

14.461 

0.767 

8.0 

10.392 

1.620 
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Table  8 

Dimensionless  stress  intensity  factor  results 
for  the  sinusoidally  varying  uniform  heat  source, 
Edge  crack  geometry 
To  =  hoL^ 


(0* 

K 

<D 

(x  10-3) 

(radians) 

0 

7.072 

-3.142 

0.5 

5.978 

-2.448 

1.0 

4.339 

-1.984 

1.5 

3.618 

-1.677 

2.0 

2.530 

-1.447 

3.0 

1.782 

-1.092 

4.0 

1.336 

-0.807 

8.0 

0.413 

0.041 
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Table  9 

Dimensionless  stress  intensity  factor  and  temperature  field  results,  o/W  =  1/6 
for  the  sinusoidally  varying  uniform  heat  source, 

Center  crack  geometry 
Xo  =  hoL2 


CO* 

K 

O 

Ti 

o 

(X  10-3) 

(radians) 

(radians) 

0 

39.801 

3.142 

0.500 

0 

1.0 

24.326 

4.179 

0.307 

0.940 

2.0 

13.962 

4.596 

0.179 

1.263 

4.0 

6.607 

-1.287 

0.0904 

1.486 

8.0 

2.467 

-0.913 

0.0419 

1.592 
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Table  10 

Dimensionless  stress  intensity  factor  and  temperature  field  results,  L/W  =  4.0 
for  the  sinusoidally  varying  uniform  heat  source 
Center  crack  geometry 
to  =  hoL^ 


to* 

K 

<D 

Ti 

(x  10-3) 

(radians) 

(radians) 

0 

9.459 

3.142 

0.500 

0.0 

0.5 

7.956 

3.834 

0.421 

0.585 

1.0 

5.775 

4.297 

0.307 

0.940 

2.0 

3.303 

4.592 

0.179 

1.263 

4.0 

1.541 

-0.814 

0.0904 

1.486 

8.0 

0.548 

0.0282 

0.0419 

1.592 
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Table  11 

Dimensionless  stress  intensity  factor  and  temperature  field  results 
for  the  antisymmetric  problem  of  sinusoidally  varying  end  temperature, 

Center  crack  geometry 


0)* 

K 

O 

Ti 

«1> 

(x  10-3) 

(radians) 

(radians) 

0 

34.319 

0.0 

0.146 

0.0 

0.5 

34.364 

0.409 

0.144 

0.287 

1.0 

34.484 

0.798 

0.137 

0.564 

1.5 

34.608 

1.151 

0.128 

0.826 

2.0 

34.674 

1.462 

0.118 

1.068 

3.0 

34.492 

1.967 

0.0972 

1.497 

4.0 

33.901 

2.349 

0.0793 

1.864 

8.0 

27.141 

3.209 

0.0373 

2.994 
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Table  12 

Dimensionless  stress  intensity  factor  and  temperature  field  results 
for  the  antisymmetric  problem  of  sinusoidally  varying  linear  heat  source. 

Center  crack  geometry. 

To  ~  hoL^ 


CO* 

K 

<D 

Ti 

(x  10-3) 

(radians) 

(radians) 

0 

4.117 

0.0 

0.0122 

0.0 

0.5 

4.054 

0.211 

0.0120 

0.208 

1.0 

3.878 

0.411 

0.0115 

0.407 

1.5 

3.629 

0.599 

0.0107 

0.589 

2.0 

3.335 

0.762 

0.00986 

0.753 

3.0 

2.757 

1.040 

0.00815 

1.024 

4.0 

2.263 

1.253 

0.00669 

1.233 

8.0 


1.111 


1.766 


0.00327 


1.718 


IV.  DIRECT  RESISTANCE  HEATING  OF  A  CENTER-CRACKED  PLATE 


An  alternative  method  for  heating  fatigue  specimens  is  to  apply  an  electric  current 
uniformly  across  the  end  of  the  specimens.  The  major  difference  between  the  direct  current 
and  heat  lamp  methods  is  that  the  insulated  crack  faces  permit  no  current  flow  and  so  the 
heat  source  generated  by  the  current  is  singular  at  the  crack  tip  although  the  temperature 
remains  bounded  there.  This  singularity  may  induce  a  different  stress  intensity  behavior 
than  the  regular  heat  source  induced  by  conventional  lamp  heating. 


Defining  the  Heat  Source 

The  heat  source  induced  by  an  electric  current  flowing  in  a  plate  is  calculated  from 
the  electric  current  density  vector,  I,  by  the  relation 


h  (x,y)  =  -^  ( I  •  I ) 


(54) 


where  p  is  the  electric  resistivity  and  k  is  the  thermal  conductivity  of  the  material.  Thus,  in 
order  to  evaluate  the  heat  source  magnitude  for  the  center-cracked  plate  problem,  we  must 
first  determine  the  resultant  electric  curent  vector  in  a  plate  of  width  W  with  a  center  crack 
exposed  to  a  far-field  electric  potential.  We  approach  this  problem  by  taking  advantage  of 
the  analogy  of  fluid  flow  around  an  infinite  cascade  of  flat  plates  and  consider  the  solution 
in  the  rectangular  region  under  investigation.  A  full  description  of  the  derivation  of  the 
electric  potential  heat  source  is  given  in  Appendix  A.  If  the  electric  potential  in  the  plate  is 
<j>,  it  is  related  to  the  electric  current  by  the  vector  expression. 


I  =  -  V  ((>/p 


(55) 


The  electric  potential  in  the  region  R  exposed  to  a  potential  drop  per  unit  length  of  <)>o  is 
derived  in  Appendix  A  as: 


<l>  =  <>o  In  [  ^  ±  2y[%  ( a  cos  y  -  P  sin  y )] 


(56) 
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where 


a  =  cos 


cosh 


2W 

B  =  sin  sinh 

2W  2W 


Co  =  cos 


71(3 

2W 


T1  = 


arctan 


-2ap 
o?-  -  al 


I 


(57) 


The  heat  source  is  proportional  to  the  square  of  the  gradient  of  <|>  and  behaves 
assymptotically  as  1/r  near  the  crack  tip  where  r  is  the  distance  from  the  crack  dp.  At  points 
distant  from  the  crack,  the  heat  source  magnitude  approaches  a  constant  value  of  <t>o- 

The  temperature  fields  induced  by  this  heat  source  are  calculated  with  the  finite 
difference  routine  described  in  detail  in  the  previous  chapter.  However,  the  singularity  of 
the  heat  source  at  the  crack  tip  may  cause  some  numerical  complications  which  were  not 
addressed  by  the  verifications  of  the  finite  difference  method  already  described.  Because  of 
this,  it  is  important  to  evaluate  the  convergence  of  the  temperature  solution  with  this  of 

behavior  in  the  forcing  function.  To  do  this,  we  consider  a  problem  in  which  the  heat 
source  term  varies  as  l/r  and  a  closed-form  solution  is  obtainable,  and  compare  the 
solution  with  the  numerically  determined  values  in  which  the  appropriate  boundary 
conditions  are  input.  The  problem  we  consider  is  the  temperature  distribution  in  a  cracked 
disk,  the  region  of  interest  readily  defined  in  cylindrical  coordinates: 

R  =(r,0lo<r<l,-7t<e<ic} 

with  the  crack  located  along  the  line  0  =  7t ,  0  <  r  <  1. 

We  seek  then,  the  solution  to  the  heat  conduction  equation  with  a  sinsuoidally 
varying  heat  source  which  behaves  as  l/r. 

V^T  =  ^  ■ 

dt  r 


(59) 
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where  the  temperature,  T  =  T(  r,  0,  t ).  We  constrain  the  temperature  on  the  disk  boundary 
to  be  zero  and  the  crack  faces  to  be  insulated.  The  solution  is  obtained  by  separation  of 
variables  technique  and  is  a  Bessel  function  of  the  first  kind,  of  zero  order,  in  r.  The 
steady  state  distribution  is  found  to  be: 


^  Xn  +  01^ 


(60) 


where  the  coefficients,  Hn,  are  defined  by  Fourier  series  representation. 


H„  = 


(61) 


the  eigenvalues,  Xn  ,  are  the  zeros  of  the  equation  Jo(Xn)  *  0,  and  Ji  is  the  first  order 
Bessel  function  of  the  first  kind. 

To  obtain  this  temperature  distribution  numerically,  we  specialize  now  to  the 
geometry  of  the  center-cracked  plate  and  input  the  appropriate  gradients  of  the  above 
distribution  along  the  vertical  edges  and  the  appropriate  temperature  values  along  the  top 
edge  as  boundary  conditions  in  the  finite  difference  code.  We  then  check  spatial 
convergence  by  performing  the  finite  difference  calculations  using  different  mesh  sizes. 
Four  mesh  sizes  of  uniform  refinement  were  employed;  1/12  (coarse),  1/24  (medium), 
1/48  (fine),  1/96  (superfine). 

Since  the  crack  tip  is  the  singular  point  of  interest,  we  check  convergence  of  the 
solution  at  that  location  and  perform  the  analysis  for  three  different  frequencies:  o  =  0, 4, 
8.  The  results  of  the  convergence  check  are  listed  in  Table  13.  The  values  are  normalized 
by  the  exact  value  from  the  closed  form  solution  for  the  cosine  temperature  component,  Tc 
and  for  the  sine  temperature  component,  Ts.  The  results  show  good  convergence  with 
errors  less  than  1  per  cent  for  the  cosine  component  and  on  the  order  of  2  per  cent  for  the 
sine  component  for  the  superfine  mesh. 
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Stress  Intensity  Results 

With  the  verification  of  the  finite  difference  method,  we  can  perform  stress  intensity 
calculations  with  the  electric  potential  heat  source  as  the  forcing  function  in  the  conduction 
equation  and  the  temperature  constrained  to  be  zero  on  the  ends.  The  results  of  these 
calculations  are  listed  in  Table  14  for  the  center  cracked  problem  for  the  normalized  crack 
length  of  1/3.  In  order  to  compare  the  stress  intensity  values,  we  again  normalize  by  the 
crack  tip  temperature.  These  data  are  plotted  in  Figure  9  and  show  that  when  a  specific 
temperature  is  maintained  at  the  center  of  the  specimen,  the  stress  intensity  factor  increases 
with  the  frequency.  This  contrasts  with  the  results  of  the  uniform  heat  problem  in  which 
the  stress  intensity  decreased  when  a  specific  temperature  is  maintained  at  the  center  line 
(Figure  5). 

The  reason  for  this  increase  is  apparent  from  the  temperature  fields  generated  by  the 
electric  current  The  electric  potential  method  of  heating  gives  a  temperature  field  which  is 
spatially  dependent  in  both  x-  and  y-  directions.  This  is  different  from  the  uniform  heating 
method  which  gives  a  one-dimensional  temperature  field  (for  the  symmetric  case).  The 
gradients  in  the  temperature  in  the  x-direcdon  increase  with  the  frequency  as  shown  in 
Figure  10  and  the  magnitude  of  the  stress  intensity  is  governed  by  the  size  of  the 
temperature  gradient 

In  order  to  compare  the  uniform  heat  stress  intensides  with  resistance  heat  stress 
intensides  for  the  case  of  a  specified  end  temperature  maintained,  we  use  the  superposidon 
procedure  described  in  Equations  (41)  through  (44)  and  the  data  of  Table  14  and  Table  3. 
The  two  examples  considered  in  Secdon  n  for  uniform  headng  are  again  invesdgated  here. 
These  are  the  temperature  at  the  crack  dp  (i)  in  phase  with  the  end  temperature  and  (ii)  out 
of  phase  with  the  end  temperature.  The  data  of  Table  3  are  used  in  conjuncdon  with  those 
of  Table  14  in  the  superposition  and  the  results  are  plotted  in  Figures  11  and  12. 
Comparing  these  results  with  the  similar  results  obtained  in  the  uniform  heat  problem,  we 
observe  that  there  is  a  difference.  First,  for  the  static  case  where  the  temperature  at  the 
center  is  in  phase  with  the  end  temperature,  the  uniform  heat  solution  gives  zero  stress 
intensity.  This  is  not  the  case  when  the  specimen  is  heated  electrically  because  of  the 
presence  of  non-^ero  temperature  gradients  in  the  x-direction.  Secondly,  the  stress 
intensity  of  electrically  heated  specimens  is  observed  to  increase  with  the  frequency,  rather 
than  decrease  as  in  the  uniformly  heated  specimens. 

Unlike  the  uniform  heat  problem,  with  the  electric  resistance  method  of  heating  we 
have  an  opportunity  to  avoid  direct  measurement  of  the  temperatures  by  analytically 
estimating  the  temperature  values  at  the  specimen  end.  For  example,  we  can  estimate  this 
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temperature  as  the  far-field  temperature  in  an  equivalently  heated  long  strip*.  The  far-field 
solution  in  this  case  is; 

T(x,y,t)=ife'§  sintot 

K  this  temperature  is  taken  as  the  end  temperature,  the  stress  intensity  factor  in  a  cracked 
plate  with  electric  current  heating  can  be  estimated  without  directly  measuring  the  plate 
temperatures^.  We  again  use  the  superposition  equations  (41)  through  (44).  The  results  of 
this  superposition,  where  the  stress  intensities  are  normalized  by  the  temperature  at  the 
crack  tip,  are  plotted  in  Figure  13  and  show  the  stress  intensity  assymptotically  approaches 
zero  as  the  frequency  approaches  zero  and  peaks  at  a  normalized  frequency  of 
approximately  4.0. 

It  is  apparent  from  Figure  13  that  the  oscillating  stress  intensity  variation  with 
frequency  behaves  differently  from  the  observed  monotonic  variation  of  the  case  when  a 
zero  end  temperature  is  prescribed  (as  in  Figure  5).  We  can  attribute  this  to  the  inverse 
dependence  of  the  magnitude  of  the  far-field  temperature  on  frequency  and  to  the 
importance  of  the  phase  lag  of  the  end  temperature  with  respect  to  the  crack  tip  temperature. 
These  two  effects  are  manifested  in  three  regions  of  the  frequency,  as  shown  in  Figure  13. 
At  very  low  frequencies,  the  end  temperature  is  quite  high  (from  Equation  (62))  and  the 
resulting  superposition  is  similar  to  the  case  of  a  harmonic  temperature  applied  at  grips  with 
no  heat  source.  Since  the  thermal  gradients  are  small  with  this  mode  of  heating,  the  stress 
intensity  factors  are  also  small. 

On  the  other  hand,  as  the  frequency  gets  large,  it  is  clear  from  Equation  (62)  that  the 
end  temperature  gets  quite  small.  In  this  case,  the  centerline  temperature  is  driven  solely  by 
the  electric  resistance  heat  source  and  the  solution  approaches  that  of  the  resistance  heat 
source  problem  with  zero  end  temperature  (evidenced  by  the  closeness  of  the  solutions  at 
frequencies  greater  than  16).  In  the  middle  of  the  frequency  regime,  the  phase  difference 
between  the  end  and  center  temperatures  plays  an  important  role.  When  the  temperatures 
are  out  of  phase,  large  thermal  gradients  exist  in  the  specimen  causing  the  stress  intensities 
to  be  quite  large. 


*This  procedure  is  less  accurate  but  has  the  advantage  that  it  requires  less  information  than  the  direct 
measurement  method.  In  other  words,  it  is  only  necessary  to  know  the  temperature  range  over  which  the 
specimen  is  cycled  in  order  to  estimate  the  stress  intensity. 

^In  practice,  the  end  grips  are  cooled  and  consequently  the  temperatures  may  be  significantly  different  from 
those  used  in  this  estimate. 
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Crack  length  scaling 

We  again  investigate  how  the  stress  intensity  scales  with  crack  length  by 
performing  the  above  calculations  with  a  normalized  crack  length  of  1/6.  The  results  of  the 
electric  potential  stress  intensity  and  crack  tip  temperature  calculations  for  the  case  in  which 
the  end  temperature  is  zero  are  listed  in  Table  15.  The  results  of  the  superposition  with  the 
far-field  temperature  data  are  shown  in  Figure  14.  Comparing  these  results  with  those  in 
Figure  11  it  is  evident  that  the  stress  intensity  again  scales  approximately  as  fll/2  at  low 
frequencies  but  the  scaling  is  frequency-dependent.  At  the  higher  frequencies  (<  8.0),  the 
stress  intensity  scales  more  closely  to  a  factor  of  a . 

Summary  of  Results 

In  this  section,  we  evaluated  the  difference  between  heat  lamp  and  electric  potential 
methods  of  cycling  the  temperature  in  a  specimen.  In  general,  the  stress  intensity  of 
electrically  heated  specimens  increases  with  frequency  rather  than  decreases  as  in  heat  lamp 
heated  specimens.  Further,  the  scaling  of  stress  intensity  with  crack  length  varies  with  the 
frequency,  scaling  as  a^f^  at  low  frequencies  and  as  a  at  higher  frequencies. 


Convergence  check  of  finite  difference  solution  for  singular  heat  source. 

Precise  value  is  1.0 


mesh  size 

o 

II 

3 

II 

3 

©=  8 

Tc 

Tc 

Ts 

Tc 

Ts 

coarse 

0.9557 

0.9374 

1.0273 

0.8913 

1.0224 

medium 

0.9720 

0.9628 

1.0046 

0.9400 

0.9998 

fine 

0.9814 

0.9850 

0.9904 

0.9803 

0.9862 

superfine 

0.9989 

1.0017 

0.9830 

1.0088 

0.9793 
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Table  14 

Dimensionless  stress  intensity  factor  and  temperature  field  results 

Electric  potential  heat  source,  a  /W=  1/3 
CT  2  2 

to  =  k  <>0  ^ 

Center  Crack 


to* 

K 

O 

Ti 

C>i 

(x  10-3) 

(radians) 

(radians) 

0 

34.560 

-3.141 

0.537 

0.0 

0.01 

34.639 

-3.126 

0.537 

0.0129 

0.1 

34.365 

-2.987 

0.533 

0.129 

0.5 

29.207 

-2.436 

0.453 

0.577 

1.0 

21.536 

-1.958 

0.330 

0.925 

1.5 

16.381 

-1.639 

0.247 

1.117 

2.0 

13.135 

-1.400 

0.193 

1.232 

3.0 

9.513 

-1.041 

0.132 

1.361 

4.0 

7.614 

-0.771 

0.098 

1.427 

8.0 

4.685 

-0.081 

0.046 

1.483 

16.0 

2.964 

0.527 

0.0225 

1.382 

32.0 

1.778 

0.980 

0.0126 

1.224 
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Table  15 

Dimensionless  stress  intensity  factor  and  temperature  field  results 

Electric  potendal  heat  source,  a  /W=  1/6 
^  ^2.2 
“Co  =  F  <l>o  L 


(D* 

K 

O 

Ti 

<I>1 

(X  10-3) 

(radians) 

(radians) 

0 

38.479 

-3.141 

0.507 

0.0 

0.01 

38.478 

-3.127 

0.507 

0.013 

0.1 

38.171 

-2.995 

0.503 

0.130 

0 

32.435 

-2.476 

0.427 

0.584 

1.0 

23.695 

-2.040 

0.311 

0.939 

2.0 

13.899 

-1.558 

0.182 

1.262 

4.0 

7.146 

-1.034 

0.092 

1.485 

8.0 

3.500 

-0.434 

0.042 

1.591 

16.0 

1.840 

0.120 

0.020 

1.574 

32.0 


1.137 


0.505 


0.010 


1.522 
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T,/To 


Figure  11.  Amplitude  of  the  stress  intensity  as  a  function  of  crack  dp  temperature  and 
dimensionless  frequency.  The  temperature  at  the  crack  rip,  Ti,  is  maintained 
by  electric  protenrial  hearing  and  is  in  phase  with  the  end  temperanire,  To- 
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Tj/To 


Figure  12.  Amplitude  of  the  stress  intensity  as  a  function  of  crack  tip  temperature  and 
dimensionless  frequency.  The  temperature  at  the  crack  tip,  Ti,  is  maintained 
by  electric  protential  heating  and  is  180®  out  of  phase  with  the  end  temperature. 

To. 
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Figure  13.  Stress  intensity  nonnalized  by  crack  tip  temperature  for  a  /W=  1/3. 

The  temperature  at  the  end  is  taken  as  the  far-field  solution  of  Equation 
(62). 


V.  IMPLICATIONS  FOR  FATIGUE  TESTING 


An  accurate  and  computationally  efficient  method  for  determining  the  effect  of 
cyclic  thermal  loading  on  crack  tip  stress  fields  has  been  developed  and  has  been  used  to 
estimate  the  effect  of  simple  thermal  transients  on  crack  tip  stress  states.  The  method  is 
generally  applicable  since  it  enables  calculation  of  stress  intensity  values  for  any  fiequency 
of  the  thermal  loading,  can  be  used  for  both  singly-  and  multiply-connected  regions,  and 
gives  a  dimensionless  result  that  applies  to  all  linear  elastic  homogeneous  materials. 

The  results  of  the  foregoing  analyses  indicate  that  stress  intensity  factors  of  cracked 
components  exposed  to  thermal  fatigue  loading  conditions  have  a  significant  dependence  on 
the  frequency  of  the  thermal  cycle.  The  dependence  is  greatest  for  electric  potential  heating 
methods  due  to  the  large  gradients  in  temperature  induced  by  the  singular  heat  source.  The 
fundamental  difference  in  the  dependence  is  that  the  stress  intensity  induced  by  uniform 
heating  decreases  with  frequency  while  the  stress  intensity  induced  by  electric  heating 
increases  with  frequency.  In  addition,  the  scaling  of  the  stress  intensity  with  crack  length 
for  electric  heating  is  found  to  depend,  to  some  extent,  on  frequency.  This  contrasts  with 
uniform  heating  in  which  the  stress  intensity  invariably  scales  as 

A  primary  goal  of  our  analysis  is  to  provide  estimates  of  the  thermal  contribution  to 
the  stress  intensity  factor  in  thermo-mechanical  fatigue  testing.  With  these  estimates,  we 
can  compare  thermally  induced  stress  intensities  to  those  induced  by  mechanical  loading. 
An  introduction  to  these  calculations  is  accomplished  by  considering  the  example  case 
described  in  Wilson  and  Warren  [21].  Samples  of  IN  100  were  thermo-mechanically 
fatigued  within  the  temperature  range  of  SOO^F  -  lOOO^F  at  thermal  frequencies  of  0.5 
cycles  per  minute,  (with  the  specimen  dimensions  and  material  data  given,  this  corresponds 
approximately  to  a  normalized  frequency,  (a*  =  2.0).  The  maximum  stress  intensity  occurs 
at  the  critical  crack  length,  roughly  one-third  the  sample  width,  and  has  a  value  of 
approximately  80  ksiV(in.).  Comparing  these  data  to  the  data  presented  in  Section  11 
(uniform  heat  problem),  a  normalized  value  of  stress  intensity  20%  that  of  the  maximum 
mechanically  induced  stress  intensity  value  reported  would  correspond  to  about  0.15. 
From  Figure  6,  if  the  centerline  temperature  is  in  phase  with  the  far-field  temperature,  the 
temperature  at  the  crack  tip  would  have  to  be  at  least  twice  as  great  as  the  magnitude  of  the 
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temperature  cycle.  If  the  two  temperatures  were  out  of  phase,  however,  the  centerline 
temperature  need  only  be  equal  to  the  far-field  temperature  in  order  to  result  in  this 
magnitude  of  stress  intensity. 

By  way  of  a  second  example,  we  consider  actual  thermal  fatigue  data  and 
superposition  to  estimate  the  stress  intensity  factor  for  another  case.  Practically,  the 
information  required  for  the  superposition  and  stress  intensity  calculations  can  be  obtained 
by  appropriately  placing  thermocouples  to  measure  the  temperature  and  phase  relations  at 
two  locations  (at  the  crack  tip  and  at  the  specimen  grips).  We  have  obtained  this 
information  for  titanium  sheet  samples  from  the  Materials  laboratories  at  the  Wright- 
Patterson  Air  Force  Base.  The  experiment  consisted  of  1.5  by  3  inch  center-cracked 
specimens  (Ti-6  Al-4  V)  which  were  heated  by  four  symmetrically  placed  lamps  at  two 
frequencies:  2.8  min  per  cycle  (21.4  cycles  per  hour)  and  7.4  min  per  cycle  (8.1  cycles  per 
hour).  Thermocouples  placed  at  both  ends  of  the  specimen  confirmed  the  symmetry  of  the 
temperature  distribution.  The  temperature  ranges  for  the  two  cases  at  the  two  locations  are 


as  follows: 

Frequency 

Crack  tip 

Grip  location 

(cy/hr) 

Tmin 

Tmax 

Tjnin 

Tmax 

8.1 

290OF 

870OF 

HOOF 

230OF 

21.4 

360OF 

870OF 

1350F 

2250F 

In  each  case,  no  measurable  phase  lag  was  noted.  Retesting  confirmed  consistency  of  the 
results. 

We  apply  this  information  and  the  data  in  Tables  3,  4  and  14  to  the  superposition 
Equations  (41)  through  (44)  to  evaluate  the  Mode  I  thermally  induced  stress  intensity  factor 
for  these  particular  cases.  For  example,  for  the  frequency  of  8.1  cy/hr,  Tq  =  tiO^F  and 
Ti  =  290OF. 

With  the  materials  data  for  Ti-6  Al-4  V : 

Elastic  modulus  (E)  16.6  xlO®  psi 

Density  (p)  0. 1 6  Ibfin^ 

Coefficient  of  thermal  expansion  (a)  5. 1  xlO'^  in^n/oF  (at  6()0OF) 

Thermal  conductivity  (k)  72.0  Btu  in/ft^/hr/oF  (at  600OF) 

Heat  capacity  (Cp)  0. 1 25  Btu/lb/®F 
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actual  estimates of  the  stress  intensities  can  be  made.  The  measured  frequencies 
correspond  to  dimensionless  frequencies  of  approximately  l.S  and  4.0.  In  order  to  make  a 
direct  comparison  between  uniform  heating  and  electric  resistance  heating,  we  use  the 
corresponding  information  for  these  frequencies  from  Tables  3,  4  and  14  to  give  the 
following  estimates  of  the  thermal  stress  intensities  induced  by  the  two  different  modes  of 
heating: 


Frequency _ Uniform  heating _ Electric  heating 

8.1cy/hr  1.76  ksiV(in.)  1.66  ksiV(in.) 

21.4cy/hr  1.52  ksiV(in.)  1.93  ksiV(in.) 

These  values  may  be  compared  to  a  fracture  toughness  of  approximately  30 
ksiV(in.)  for  this  particular  alloy.  The  results  above  are  consistent  with  the  previous 
observation  that  the  dependence  of  stress  intensity  factor  on  frequency  is  different  for  the 
two  modes  of  heating  and  indicate  that  for  higher  frequencies,  much  larger  errors  in  the 
overall  stress  intensity  calculations  may  actually  be  occurring  when  electric  resistance 
heating  is  used  to  cycle  the  specimen  temperature  in  thermo-mechanical  fatigue. 

It  is  noted  that  the  error  from  the  dynamic  temperature  distribution  contributes  only 
a  part  of  the  total  error  stemming  from  the  incorporation  of  thermal  cycling  in  thermo¬ 
mechanical  fatigue.  This  is  because  the  temperature  at  the  center  does  not  oscillate  about 
the  same  temperature  as  the  end  temperature.  In  fact,  a  static  temperature  gradient,  which 
we  approximate  as  parabolic  in  the  y-direction,  exists  in  addition  to  the  dynamic  oscillations 
about  this  gradient.  Since  the  center  of  the  specimen  is  measured  as  hotter  than  the  end  of 
the  specimen,  the  stress  intensity  contribution  from  this  static  thermal  gradient  is  negative. 
We  can  estimate  this  contribution  using  the  static  data  of  Table  4  and  the  materials 
information  cited  above  as  12.5  per  cent  and  15  per  cent  of  the  fracture  toughness  for  the 
dimensionless  frequencies  of  4  and  1.5,  respectively. 

Because  the  errors  stemming  from  the  dynamic  part  of  the  temperature  distribution 
oscillate,  they  can  either  amplify  or  attentuate  the  total  error  induced  by  thermal  cycling.  In 
other  words,  since  the  static  temperature  yields  a  negative  stress  intensity  and  the  error  in 
the  dynamic  pan  can  be  positive  or  negative,  summing  the  negative  errors  could  give  an 
overall  error  as  much  as  about  25  per  cent  of  the  fracture  toughness  valm.  This  means  that 


stress  intensity  values  calculated  are  estimates  in  light  of  the  assumption  that  the  higher  order  terms 
•n  the  Taylor's  series  expansion  are  considered  small  in  comparison  with  the  fust  term  and  the  aitalysis  was 
performed  in  plane  strain  while  plane  stress  samples  were  us^  fcv  testing. 
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the  load  cycle  could  be  significantly  lower  than  what  is  actually  recorded.  Further,  since 
the  applied  maximum  cyclic  load  often  reaches  only  one-third  to  one-half  the  fracture 
toughness  value,  the  error  in  the  loading  stemming  from  the  cyclic  thermal  loads  could,  in 
fact,  be  significantly  greater  than  25  per  cent 

As  discussed  in  Section  IV,  the  electric  resistance  method  of  heating  gives  an 
opportunity  to  avoid  direct  measurement  of  the  temperatures  by  analytically  estimating  the 
temperature  values  at  the  specimen  end  as  the  far-field  temperature  in  an  equivalently  heated 
long  strip.  Equation  (62).  If  this  temperature  is  taken  as  the  end  temperature,  the  stress 
intensity  factor  in  a  cracked  plate  with  electric  current  heating  can  be  estimated  without 
directly  measuring  the  plate  temperatures. 

It  is  interesting  to  compare  this  method  of  estimating  thermally  induced  stress 
intensity  factors  with  the  method  provided  in  [35],  in  which  it  is  proposed  that  a  static 
model  problem  be  used  to  estimate  the  stress  intensities  induced  by  direct  resistance 
heating.  The  results  of  this  earlier  analysis  yielded  axisymmetric  temperature  variations  in 
the  vicinity  of  the  crack  tip  as  T(r)  =  Tb+  Tc(l-r)  which  gives  a  resulting  Mode  I  stress 
intensity  factor  as 

K,  =  0.0914 
I  l-v  ^ 

where  Tb  is  a  constant  temperature  along  a  circular  boundary  at  a  radius  of  the  half  crack 
length,  a,  r  is  the  dimensionless  distance  from  the  crack  tip  (normalized  to  a )  and  Te  is  a 
temperature  which  is  a  function  of  the  electrical  and  thermal  conductivities  of  the  specimen, 
the  applied  voltage  drop  per  unit  length  and  the  crack  length: 

With  the  test  data  and  INIOO  materials  properties  information  presented  in  [35],  the  Mode  I 
stress  intensity  factor  at  the  highest  voltage  is  calculated  to  be  3.8  ksiV(in.). 

We  can  estimate  the  Mode  I  stress  intensity  with  the  far-field  approximation  of  the 
end  temperature  described  above.  If  a  complete  thermal  cycle  is  achieved  in  one  minute 
then  the  normalized  frequency  corresponding  to  this  rate  of  heating  is  approximately  4. 
Using  the  data  of  Table  14  and  the  superposition  procedure  described  in  Section  II,  Kj  is 
calculated  to  be  3.8  ksi‘\/(in.).  Although  a  comparison  of  the  model  problem  estimate  with 
this  result  shows  a  very  close  correlation  for  this  particular  case,  the  results  depicted  in 
Figure  13  show,  in  general,  the  model  problem  will  give  results  within  an  order  of 
magnitude  of  that  predicted  with  the  dynamic  analysis  discussed  here. 

It  has  been  shown  that,  with  the  results  presented,  estimates  of  the  thermal 
contribution  to  stress  intensity  in  thermo-mechanical  fatigue  testing  can  be  derived  for  a 
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wide  variety  of  symmetric  loading  condidons.  However,  these  calculadons  are  not  limited 
to  symmetric  temperature  distributions.  With  the  appropriate  superposition  and 
measurements  of  temperature  at  several  locations  within  the  sample,  we  can  derive  both 
Mode  I  and  Mode  n  stress  intensity  values  for  asymmetric  thermal  loadings.  It  is  of 
interest  to  evaluate  how  large  the  antisymmetric  part  of  the  temperature  field  has  to  be  in 
Older  for  the  Mode  n  stress  intensity  to  be  significant  compared  to  the  Mode  I  value.  In  an 
example  situation  of  a  normalized  frequency  of  2.0,  in  order  for  the  Mode  I  and  Mode  II 
stress  intensities  to  be  comparable,  the  antisymmetric  end  temperature  would  have  to  be 
five  times  greater  than  the  symmetric  end  temperature  and  the  center  temperature  five  times 
less. 

Geometrical  variations  have  been  investigated.  We  have  found  that  edge-cracked 
specimens  result  in  much  lower  stress  intensity  values  (by  a  factor  of  approximately  five) 
than  the  center-cracked  specimens.  This  is  probably  due  to  the  lack  of  restraint  in  bending 
that  exists  in  the  edge-cracked  geometry  that  does  not  exist  in  the  center-cracked  plate. 

In  conclusion,  the  results  of  this  analysis  indicate  that  significant  stress  intensities 
may  be  induced  by  thermal  fatigue  loading.  Further,  the  stress  intensity  values  depend  on 
the  mode  of  heating  used  to  cycle  the  specimen.  Samples  heated  by  direct  resistance 
heating  are  subjected  to  a  singular  heat  source  at  the  crack  tip  which  results  in  the  stress 
intensity  factor  increasing  with  frequency.  This  behavior  is  opposite  that  of  uniformly 
heated  specimen  and  indicates  that  care  should  be  taken  when  resistance  heating  is  used  at 
relatively  high  frequencies.  The  results  obtained  in  the  preceding  sections  may  be  used  in 
conjunction  with  experimental  data  to  improve  predictions  of  stress  intensity  factors  in 
laboratory  fatigue  tests  and  in  actual  fatigue  service  environments. 


VI.  A  PILOT  STUDY  OF  COMPOSITE  MATERIALS 


In  the  interest  of  understanding  the  effects  of  thermal  transients  in  composite 
materials,  an  appropriate  extension  of  the  path  independent  integral  approach  is  its 
application  to  bimaterial  composites  in  which  cracks  occur  normal  to  the  bimaterial 
interface.  We  select  this  particular  geometry  for  investigation  because  of  its  physical 
interest:  it  can  correspond  to  the  problem  of  a  crack  in  the  matrix  or  in  a  fiber,  turning  at  or 
penetrating  the  interface  between  the  fiber  and  the  matrix;  and  because  the  corresponding 
eigenvalues  are  real.  Real  eigenvalues  permit  the  interpretation  of  the  eigenfunctions  in 
terms  of  standard  physical  arguments  since  there  is  no  oscillating  singularity  at  the  crack 
tip. 

The  isothermal  problem  of  a  composite  with  a  crack  normal  to  the  interface  has  been 
studied  by  Hilton  and  Sih  [28]  for  the  case  in  which  the  crack  does  not  extend  to  the 
interface,  and  by  Cook  and  Erdogan  [29]  and  Gupta  [30]  for  the  limiting  case  where  the 
crack  ends  at  the  interface.  In  each  case,  the  method  of  solution  uses  integral  transform 
techniques  to  extract  the  stress  intensity  factors  for  different  material  combinations. 

In  our  approach,  we  consider  the  relative  importance  of  two  thermal  loading  effects, 
those  induced  by  dissimilar  thermal  diffusivities  and  dissimilar  thermal  expansions.  In 
performing  these  parametric  studies  we  are  able  to  use  the  same  techniques  as  those 
developed  for  the  homogeneous  case.  Finite  difference  approximations  will  be  derived  for 
the  temperature  fields,  which  now  depend  on  the  ratio  of  the  thermal  diffusivities  of  the  two 
materials.  Finite  element  approximations  will  be  obtained  for  the  stress  and  displacement 
vectors  which  now  depend  on  the  ratio  of  the  thermal  expansions  and  the  elastic  moduli  of 
the  two  materials.  Finally,  a  path  independent  integral  will  be  used  to  extract  the  stress 
intensity  factors  from  the  appropriate  stress,  displacement  and  temperature  fields. 

Formulation 

In  order  to  investigate  the  effects  of  transient  thermal  loadings  on  crack  tip  stress 
fields  in  composite  materials,  we  follow  Okajima  [31]  in  formulating  the  following  set  of 
problems.  First  we  consider  the  rectangular  region  in  the  x,  y  plane  illustrated  in  Figure 
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15.  By  virtue  of  the  symmetry  of  this  configuration,  we  are  able  to  limit  the  analysis  to  the 
region  R:  which  is  composed  of  two  subregions  R  i  and  R  2- 

R  =Rj  U/?2. 

=  ( (x,y)  I -W<x<0,0<y<L} 

J  (63) 

/?2  =  1  (x.y)  I  0<x<W,  0<y<L} 


Once  again  we  seek  the  two-dimensional  thermoelastic  stress  and  displacement 
vectors  satisfying  the  plane  uncoupled  quasi-static  equations  of  thermoelasticity. 
Reiterating  the  earlier  formulation,  we  have  the  following  equations  for  bimaterials.  First, 
the  temperature,  T(x,y,t)  is  defined  by  the  heat  conduction  equation 


T  (x,y,t)  =  ^  (x.y.t)  -  H  (x,y,t) 

I 


(64) 


on  R  i  ,  i  =  1,2,  and  boundary  conditions  which  will  be  discussed  later.  Secondly,  the 
stress  vector  satisfies  the  plane  equations  of  equilibrium  on  R  in  the  absence  of  a  body 
force  field  (Equations  (7)),  and  the  stress  and  displacement  vectors  satisfy  the  plane-strain 
stress  displacement  relations: 


(65) 


on  /?,  ,  i  =  1,2.  The  stress  and  displacement  vectors  satisfy  traction-free  boundary 
conditions  on  the  outer  edges  and  along  the  crack  face: 


=  tjy  =  0  on  X  =rtW ,  (0  <  y  <  L) 

Oy  =Txy  =  0ony=L,  (-W<x<W) 

Oy  =  Txy  =  0  on  y  =  0,  (-W  <  x<  0) 


(66) 


and  symmetry  conditions  along  the  positive  x-axis: 
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v  =  Txy=0  ony  =  0,(0<x<W) 


(67) 


Finally,  the  two  regions  are  considered  perfectly  bonded  and  so  must  obey  continuity 
requirements  at  the  interface: 


and 


lim  ^  Oj  =  lim  Ox  ,  lim  t^v  =  ^  '^xv  (0  <  y  <  L) , 

X  -*o  x-tcr  x-iO  ^  x-»(r 


(68) 


lim  u  =  lim  u  ,  lim  v  =  lim  v  (0  <  y  <  L) 
X  x-i(r  x^*  x-*(r 


(69) 


where  0+  implies  x  >  0,  and  0*  implies  x  <  0. 

Again,  we  refer  to  the  stress  intensity  factor  for  evaluating  the  importance  of 
transient  temperature  fields  on  crack  tip  stresses.  Herein  we  deviate  from  homogeneous 
behavior.  For  cracked  non-homogeneous  solids  with  different  elastic  moduli  the 
eigenvalue  is  no  longer  equal  to  1/2.  If  we  denote  X  as  the  singular  eigenvalue 
characterizing  the  only  singular  stress  field  at  the  crack  tip,  the  stress  intensity  factor  for 
composite  materials,  Kx  is  defined  as 

K.  =  lim^  x^’^Ov  on  y  =  0 

^  (70) 

Evaluation  of  the  eigenvalues  and  eigenfunctions 

In  order  to  determine  the  eigenvalues  for  the  present  case,  we  refer  to  the  analysis 
of  Dempsey  and  Sinclair  [32]  and  specialize  to  our  geometry.  Substituting  0i=  7t/2  and 
02=  n  into  Equation  l-A-4  on  page  323  of  their  analysis,  we  obtain  the  transcendental 
equation  for  X  as  a  function  of  material  constants  pi  and  Vj,  i  =  1,2: 

sin2^  (|3^-l)  +  X^P-p^-a+ap)+^|^  =0 


In  (71),  a  and  P  are  defined  as  combinations  of  Pi  and  Vj, 


(71) 
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“-^2^1+^1012  •  P-“  ^  ^2®1+^^1™2 


(72) 


In  plane  strain,  nii  =  4  (1-  Vi ).  The  composite  material  parameters  were  first  defined  in  this 
way  by  Zak  and  Williams  [33].  It  is  noted  that  Equation  (71)  is  identical  to  the 
transcendental  equation  derived  in  the  analysis  of  Cook  and  Erdogan.  The  limiting  values 
of  a  and  p  are  determined  by  the  physical  restrictions  of  Vj  (0  <  Vj  <  0.5),  and  iii  (0  <  p-i  < 
oo).  Substituting  these  values  into  (72)  gives  -1  <  a  <  1,  and  (a-1)  <  4p  <  (a+1). 
Several  eigenvalues  satisfying  (71)  for  this  range  of  material  parameters  are  shown  in 
Figure  16. 

To  determine  the  eigenfunctions  associated  with  the  eigenvalues,  we  turn  to  the 
work  of  Williams  [4].  The  results,  best  presented  in  the  cylindrical  coordinate  system  ( x  = 
r  cos  0,  y  =  r  sin  0 ),  are  reported  here. 


oj  = - - [  f  ••  (0;X)  +  (X+1)  f.  (0;>.)] 

/2k  X  (X+1)  b 

oe'  =  (>.+1) f, (e;>.)] 

/2k  X  (X+1)  b 

K,  r^~‘ 

tre*  =  -  [-Xf.'(0;X)] 

/2k  X  (X+1)  b 


^  A - [ .  (X+1)  f  (0;X)  +  ^  m.  g.’  (0;X)  ] 

2/2jcX(X+l)biij  ^ 

ue ‘  -  [  -f ’(e;X)  +  m  (X-1)  g.  (0;X)] 

2/27cX(X+l)bPj  ‘  ^  ‘ 


a3) 


in  i ,  i  =  1,2.  The  prime  denotes  differentiation  with  respect  to  0  and  the  functions, 
fi(0;X)  and  gi(0;X)  are 

f  (0;X)  =  bjjSin  (X+1)0  +  b^jCos  (X+1)0  +  bj^sin  (X-1)0  +  b^^cos  (X-1)0 

g.  (0;X)  =  7^  I  -^i3  (^-1)9  +  ^i4  sin  (X-1)0] 

*  1 


05) 
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where  b  =  b22+b24.  and  the  constants  by  ,  i=l,2,  j  =  1,2, 3,4,  are  determined  from  the 
traction  and  displacement  boundary  conditions: 


^21  ~  ^23  “  ®  (from  symmetry) 

^24=  1 

b22  =  B 

bj4  =  ^  {-  (X  +  1)[s2(4  -  3)->.]b22+  [s^4>ls2  -  5X  +  \)-X\  X]b^^ 


(76) 


b,,  = - - 

2sc(2Ul) 


[  2(2X+1)  s^-  2X]  +  [  2(k+\)  i-X-l]  b^ 


-(2Xs^-X+l)b24  ► 

bj2  =  -4  sc  bj3  -  (1-  4  s^)  bj4  +  2  s^  ( b^j-  b^^ ) 

^11  ““C  (^12"^14‘'’22‘'’^24)‘*'^13 
A  =  s^-X^ 

„  [(2p+2)X-2p+2Is2  +  2(a-P)xV  (2p+a-l)X^-(a+l)X 

2(j3+l)(X+l)s^  +  2(a-P)X  +(3a-2p+l)X  +  (a+l)X 

s  =  sin^^ 

2 


(77) 


Finally,  the  path  independent  integral  requires  the  divergence  of  the  displacement 
eigenvectors,  i=I,2: 

The  complementary  eigenfunctions,  a‘  *  and  u'  * ,  i  =  1,2,  arc  obtained  by  substituting  -X 
for  X  and  a  scaling  factor,  Kx*.  for  Kx,  in  Equations  (73)  through  (77). 

We  follow  the  procedure  outlined  in  Section  HI  for  defining  the  appropriate  path 
independent  integral  which  extracts  the  stress  intensity  factor.  By  virtue  of  the  continuity 
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requirements  at  the  interface,  the  path  independent  integral  is  that  of  Equation  (31)  in  which 
K  is  replaced  by  Kx.  and  the  homogeneous  complementary  eigenfunctions  are  replaced  by 
the  complementary  functions  just  developed: 

2 

Kj^  =  ^{  J[(OxU‘%x^v^*-a,‘*u-T4rv)n, 
i=l  z 

j 

+  (  Oy  V*  +  Xjy  U  -  Oy  V  -  U  )  Oy  ]  dS 

+  pJt(u‘*.*  +  v‘*.y  )  dA  ) 


where  the  integration  paths,  Ei ,  i  =  1,2  ,  refer  to  the  boundaries  not  including  the  crack 
faces,  of  the  regions  ,  i  =1,2  ,  respectively.  The  scaling  factors  ,  Kx*.  are  determined 
by  substituting  the  eigenfunctions  themselves  into  (79),  letting  T  =  0,  and  performing  the 
integration. 

To  verify  the  accuracy  of  the  established  eigenfunctions  and  the  path  independent 
integral  approach  we  analyze  a  check  problem  which  utilizes  the  most  strongly  singular 
eigenvalue  used  in  the  intended  application.  We  follow  Okajima  [31]  in  formulating  the 
check  problem  which  superimposes  eigenfunctions,  adjusting  their  participation  so  that  a 
reasonable  approximation  of  the  boundary  conditions  used  in  the  application  results.  The 
forms  of  the  superposition  are  those  of  (73)  with  Kx  being  replaced  by  the  constants  Kx 
and  Cz  corresponding  to  the  eigenvalues  and  Xz-  Specifically,  we  choose  the  material 
combination  corresponding  to  the  most  singular  eigenvalue,  the  combination  giving  the 
ratio  of  elastic  moduli  of  Ez/Ei  =  0.14.  The  first  two  eigenvalues  associated  with  this 
combination  are 

=  0.2816,  ^2  =  1.1077 


In  combining  the  first  two  eigenfunctions  to  approximate  a  uniaxial  tension  loading  on  the 
region  of  interest,  we  take  the  participation  of  the  first  eigenvalue  field ,  Kx  =  1.0  and  the 
participation  of  the  second,  Cz  =  0.25.  The  traction  conditions  along  the  edge,  then,  are 
determined  as  the  sum  of  these  two  eigenfunctions  with  the  appropriate  participation  factor. 
The  loading  conditions  for  the  check  problem  are  shown  in  Figure  17. 
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To  examine  convergence  of  the  solution  we  perform  the  analysis  on  three  uniformly 
refined  grids.  A  full  description  of  the  Hnite  element  grids  used  is  contained  in  Appendix 
B.  The  results  of  the  analysis  demonstrate  acceptable  performance  on  the  test  problem. 
The  value  of  Kx.  for  the  coarse,  medium  and  fine  grids,  are  respectively  1.061, 1.029  and 
1.018  (the  exact  value  is  1.0). 

Problem  D^nition 

The  goal  of  this  analysis  is  to  establish  the  importance  of  thermal  diffusivity 
ntismatches  relative  to  mismatches  in  thermal  expansion  rates.  In  order  to  accomplish  this, 
we  pose  two  problems.  In  the  first  we  consider  the  ratio  of  oE  in  each  material  to  remain 
unity  while  the  ratio  of  thermal  diffusivity  changes.  In  this  case,  a  frequency  dependence 
will  be  established.  The  eigenvalue  will  be  limited  to  the  case,  X  =  0.5,  since  the  the  elastic 
moduli  of  the  materials  are  equal.  The  heating  conditions  imposed  are  a  uniform  heat 
source  varying  as  hocos(cot)  with  all  edges  insulated.  Although  the  static  case  (to  =  0 ) 
gives  an  infinite  temperature  distribution,  we  will  look  at  the  behavior  of  the  stress  intensity 
factor  as  a  function  of  frequency  and  study  the  limit  as  the  frequency  approaches  zero.  In 
the  second  problem  we  consider  the  ratio  of  thermal  diffusivities  of  the  two  materials  to  be 
unity  and  evaluate  the  stress  intensity  factor  for  several  ratios  of  thermal  expansion  and 
elastic  moduli,  since  it  is  the  combination  of  these  material  constants  that  gives  rise  to  the 
magnitude  of  the  thermal  stresses  in  a  material.  In  this  problem,  no  frequency  effects  will 
result  since  the  temperature  is  uniform  throughout  the  composite.  Once  these  preliminary 
relationships  are  established,  we  will  evaluate  the  importance  of  diffusivity  mismatches  and 
expansion  ntismatches  in  some  currently  used  composites. 

Thus,  in  addition  to  the  constraints  imposed  on  the  stress  and  displacement  vectors 
presented  in  Equations  (65)  through  (71),  the  following  conditions  on  the  temperature  field 
satisfying  Equation  (63)  will  be  imposed.  First,  the  edges  of  rectangular  geometry  are 
taken  to  be  thermally  insulated. 

T,x  =0  on  x  =  0andon  x  =  W,  0<y<L 
T,y  =  0  on  y  =  0  and ony  =  L,  0<x<W 

Secondly,  either  the  heat  generation  term  behaves  harmonically  and  the  ratio  of  thermal 
diffusivity  varies  between  0.1  and  10, 


H(x,y,t)  =  ho  cos  cot,  ajEi/  a2E2  =  1, 0.1  <  Di/  D2  <  10,  Vi  =  V2.  (82) 
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or  the  temperature  remains  uniform  throughout  the  composite  and  the  ratio  of  product  of 
thermal  expansion  coefficient  times  the  elastic  modulus  varies  between  the  limits  0.1  and 
10. 

T(x,y,t)  =  To.  Di/  D2  =  1, 0.02  <  a\Ei/  a2E2  <  50,  vi  =  V2.  (83) 

The  limits  on  the  ratios  of  material  constants  have  been  chosen  to  approximate  to  first  order 
those  of  actual  composite  materials.  Regardless  of  the  material  combination,  we  have  taken 
Poisson's  ratio  to  be  the  same  in  each  material  and  equal  to  0.3.  This  assumption 
simplifies  the  calculations  required  for  determining  the  eigenvectors  since  their  dependence 
is  now  reduced  to  one  parameter;  the  ratio  of  the  elastic  moduli.  In  general,  for  the 
composite  materials  under  consideration,  Poisson’s  ratio  varies  very  little  between  the  two 
materials,  ranging  between  0.25  and  0.35. 

Thermal  Diffusivity  Effects 

Following  the  method  of  solution  outlined  in  Section  II,  Equations  (14)  through 
(16),  we  represent  the  temperature,  T(x,y,t)  as  the  sum  of  harmonic  components  in  which 
(0  is  the  frequency  of  the  thermal  excitation  and  Ts(x,y)  and  Tc(x,y)  are  the  magnitudes  of 
their  respective  components.  The  solution  to  the  heat  conduction  Equation  (63)  with  the 
imposed  conditions  of  (81)  and  (82)  for  the  limiting  case,  D2  =  Di  =  D,  is  straightforward: 

T  (x,y,t )  =  — -  -f-  sin  o)t  for  to  ^  0,  and 

W"  (84) 

D 

T(x,y,t)  =  — r-t  for  t0  =  0. 

w  (85) 

Clearly,  Ts  and  Tc  are  undefined  for  the  case  the  frequency  is  zero  and  are  spatially  constant 
for  values  of  frequency  not  equal  to  zero  ( Ts  =  D/CW^co),  Tc  =  0 ).  The  resulting  stress 
intensity  factors,  Ks  and  IQ,  for  this  limiting  case  are  zero  since  the  plate  is  not  restrained 
from  uniform  expansion. 

For  the  general  case  of  varying  D2/D1 ,  the  temperature  solution  is  more  complex. 
Finite  difference  solutions  have  been  obtained  for  values  of  D2/D1  =  0.1, 0.5, 5.0, 10.0  for 
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a  range  of  dimensionless  frequencies  from  0.03  to  2.55^  ^  These  results  are  pictured  in 
Figure  18  in  which  the  temperature  field  is  the  amplitude  of  the  harmonic  ccHnponents,  T  = 
(Ts^  +  and  has  been  normalized  by  the  temperature  amplitude  at  the  crack  tip,  Ti. 
Because  of  the  insulation  conditions  on  the  y-faces  of  the  plate,  no  temperature  variation  in 
the  y-direction  exists  and  the  solutions  illustrated  in  Figure  18  hold  for  all  values  of  y. 
From  these  figures,  a  general  trend  of  increasing  temperature  gradient  with  frequency  is 
observed  and  as  the  frequency  approaches  zero,  the  temperature  approaches  a  flat 
distribution.  The  variation  of  crack  tip  temperature  with  frequency  for  the  different  values 
of  the  diffusivity  ratio  is  shown  in  Figure  19.  Here,  the  dimensionless  crack  tip 
temperature  has  been  normalized  by  the  value  hoW^,  ho  is  the  magnitude  of  the  heat  source 
and  W  is  the  plate  width.  The  crack  tip  temperature  tends  to  increase  to  an  unknown  limit 
as  the  frequency  approaches  zero  and  its  magnitude  also  increases  with  the  diffusivity  ratio. 

From  these  temperature  distributions,  stress  intensity  factors  were  calculated  using 
the  same  procedures  detailed  in  Sections  n  and  IE.  The  finite  element  method  was  used  to 
compute  the  stress  and  displacement  fields  associated  with  the  temperature  fields  and  the 
Mode  I  path  independent  integral  (Equation  (31))  was  used  to  calculate  the  stress  intensity 
factors  from  these  fields.  Uniform  grid  refinement  was  used  for  all  calculations.  The  finite 
element  grids  are  shown  in  Appendix  B.  Convergence  behavior  of  the  results  was  good 
and  several  paths  of  integration  were  used  in  the  computations  to  demonstrate  path 
independence.  The  results  of  the  calculations  are  shown  in  Figure  20.  In  these  figures,  the 
stress  intensity  factor  represented  is  the  amplitude  of  the  harmonic  components:  K  =  (Ks^ 
+  Kc^)^^  and  has  been  nondimensionalized  by  the  crack  tip  temperature  amplitude,  Ti: 


K  = 


K, 


aET 


1-v 


IJm 


(86) 


The  trend  of  increasing  stress  intensity  with  frequency  follows  from  the  temperature 
results:  increasing  thermal  gradients  with  frequency.  It  is  also  noted  that  as  the  frequency 
approaches  zero,  so  does  the  stress  intensity. 

Thermal  Expansion  Effects 


^  ^  The  dimensionless  frequency,  (o*  ,  is  now  normalized  by  the  diffusivity  of  material  1  and  thus  is  related 
to  the  frequency,  (O,  by  the  relation  to*  =  to  L^/fDj  w). 
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To  investigate  the  effects  of  thermal  loads  in  composite  materials  with  mismatches 
in  thermal  expansion  and  elastic  modulus,  we  flrst  consider  static,  uniform  temperature 
fields.  Performing  this  preliminary  study  enables  us  to  compare  the  importance  of 
frequency  effects  relative  to  static  loading  effects  since  mismatches  in  expansion  and  elastic 
modulus  do  not  affect  the  temperature  field  of  the  composite.  Thus,  we  take  the 
temperature  in  the  composite,  T(x,y)  to  be  a  constant,  Tq,  and  the  heat  source,  H(x,y)  to  be 
zero.  This  field  satisfies  the  heat  conduction  Equation  (64)  with  the  conditions  of  thermally 
insulated  boundaries  (Equations  (81))  described  in  the  problem  definition.  Further,  we 
vary  the  ratio  of  the  product  of  themal  expansion  and  elastic  modulus  between  the  limits  of 
0.02  to  50,  as  described  in  (83). 

With  this  temperature  field,  stress  intensity  factors  were  calculated  using  the  finite 
element  method  to  compute  the  stress  and  displacement  fields  and  the  path  independent 
integral  of  Equation  (79).  The  case  of  uniform  elastic  modulus  (E2  =  Ei)  is  interesting 
because  the  corresponding  eigenvalue  is  1/2.  Results  of  the  stress  intensity  calculations  are 
shown  in  Figure  21  for  this  case.  The  results  show  that  the  stress  intensity  factor  is 
linearly  related  to  the  difference  in  the  thermal  expansion  coefficients  of  the  two  materials 
and  is  zero  when  the  expansion  coefficients  are  equal.  From  these  data,  the  stress  intensity 
factor  is  calculated  to  be 

ET 

K,  =  0.28  (a,  -  a- )  ina 

I  1-v  1  2  (87) 


For  the  case  of  equal  thermal  expansion  and  varying  elastic  modulus,  no  stress  intensity 
exists  since  this  corresponds  to  the  case  of  uniform,  unrestrained  expansion.  When  the 
expansion  coefficient  is  not  uniform,  however,  the  stress  intensity  factor  can  become  quite 
large.  Results  of  these  calculations  are  shown  Figure  22  along  with  the  eigenvalues 
associated  with  each  elastic  modulus  ratio  used  in  the  computations.  The  stress  intensity 
values  are  normalized  to  the  properties  of  material  1,  the  material  containing  the  crack,  and 
the  uniform  temperature  in  the  plate: 


K  = 


K, 


l-v 


(88) 
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From  Figure  22,  we  find  that  when  the  thermal  expansion  coefficient  is  larger  in  the 
uncracked  material  (a2/ai  >  1)  the  normalized  stress  intensity  is  greater  than  zero,  in  the 
crack-opening  mode  and  is  less  than  zero  (crack  closure)  when  the  coefficient  is  larger  in 
the  cracked  material. 

Summary  and  Discussion  of  Results 

In  this  section  we  have  extended  our  approach  of  estimating  thermal  fatigue  stress 
intensity  factors  to  a  composite  material  by  performing  a  pilot  study  of  bimaterial  composite 
with  a  crack  normal  to  the  interface.  Comparing  the  effects  of  diffusivity  and  expansion 
from  these  preliminary  studies  indicate  that  frequency  effects  are  an  order  of  magnitude 
lower  than  those  induced  by  variations  in  elastic  modulus  and  expansion.  In  order  to 
investigate  this  further,  we  have  chosen  two  composite  systems  which  are  currently  being 
considered  for  use  in  thermal  fatigue  environments.  These  are  composites  of  titanium  with 
silicon  carbide  fibers  and  glass  with  silicon  carbide  fibers.  Material  propenies  for  titanium, 
a  high  temperature  glass  trademarked  VYCOR,  and  silicon  carbide  fibers  are  listed  in  Table 
16.  We  consider  two  separate  cases  for  each  composite:  the  crack  located  in  the  matrix  (Ti 
or  glass)  and  the  crack  located  in  the  SiC  fiber.  Ratios  of  the  material  parameters  and  the 
associated  eigenvalues  for  these  four  cases  are  listed  in  Table  17.  We  have  combined  the 
thermal  diffusivity  problem  now  with  the  thermal  expansion  problem  and  have  calculated 
the  stress  intensity  factors  induced  by  cyclic  temperature  fields  where  the  heat  source  varies 
as  hgcos  (Ot  and  the  elastic  materials  properties  (E  and  a  )  are  different  for  the  two 
materials.  Results  of  these  calculations  are  shown  in  Figures  23  and  24.  The  stress 
intensity  factors  in  these  figures  have  been  normalized  by  the  amplitude  of  the  crack  tip 
temperature,  Ti ,  and  elastic  properties  of  the  material  with  the  crack,  material  1  (Equation 
(88)  with  To  replaced  by  Ti).  These  figures  iUustrate  that  for  higher  frequencies,  the  effect 
of  frequency  on  the  stress  intensity  amplitude  can  be  quite  high.  In  the  composite  of 
titanium  matrix  composite,  cyclic  temperature  fields  reduce  the  stress  intensity  regardless  of 
the  location  of  the  crack  while  different  behavior  is  observed  in  the  composite  with  the 
glass  matrix.  In  this  case,  low  frequencies  cause  the  stress  intensity  factor  to  increase 
above  the  static  value.  This  behavior  continues  at  higher  frequencies  if  the  crack  is  located 
in  the  fiber  while  the  stress  intensity  decreases  at  higher  frequencies  if  the  crack  is  located 
in  the  matrix  material. 

If  we  consider  actual  fatigue  test  frequencies  of  the  order  of  1  cycle  per  minute,  the 
stress  intensity  relative  to  the  static  value  can  be  estimated.  These  results  are  listed  in  Table 
18.  The  effect  at  this  frequency  can  be  as  great  as  increasing  the  stress  intensity  by  18  per 
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cent  (for  glass-SiC  systems  with  the  crack  located  in  the  SiC)  or  decreasing  the  stress 
intensity  by  45  per  cent  (for  Ti-SiC  systems  with  the  crack  located  in  the  Ti)  or  as  small  as 
a  2  per  cent  decrease  (for  glass-SiC  systems  with  the  crack  in  the  glass). 
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Table  16 

Material  properties  for  example  composites  12 


Material 

Heat  capacity 

(Co) 

Density 

(P) 

Thermal 

Conductivity 

(k) 

Thermal 

Expansion 

(a) 

Elastic 

Modulus 

(E) 

cal 

cal 

lO^^AC 

g°C 

g/cc 

cms®C 

in 

lO^psi 

Titanium 

0.125 

4.50 

0.0052 

8.5 

16.8 

Silicon  Carbide 

0.186 

3.22 

0.0400 

5.0 

70.0 

VYOOR 

0.180 

2.18 

0.0033 

0.75 

9.6 

l^The  diffusivity,  D,  is  related  to  the  heat  capacity,  density  and  thermal  conductivity  by  the  relation,  D  k 

/(pCp). 
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Table  17 

Ratios  of  material  properties  and  the  associated  eigenvalue  for  example  composites 


Material  1-  Material  2 

^2 

«2 

^2 

E, 

X 

Ti  -  SiC 

7.23 

0.59 

4.17 

0.6225 

SiC  -  Ti 

0.14 

1.70 

0.24 

0.3419 

Glass  -  SiC 

7.95 

6.67 

7.29 

0.6538 

SiC  -  Glass 

0.13 

0.15 

0.14 

0.2793 
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Table  18 

Stress  intensity  increase  at  a  frequency  of  1  cycle/min  for  example  composites. 


Composite 

Crack  location 

Normalized 

Frequency,  co* 

AK 

Ti  -  SiC 

Ti 

23.3 

-17% 

SiC  -  Ti 

SiC 

1.24 

-45% 

Glass  -  SiC 

Glass 

25.6 

-2% 

SiC  -  Glass 

SiC 

1.24 

+18% 

material  1 


material  2 


axis  of  symmetry 


Figure  15.  Geometry  of  the  composite  problem. 


Figure  17.  Traction  conditions  for  the  test  problem. 
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9 


0  2  4  6  8 

Frequency 


Figure  19.  Variation  of  the  crack  tip  temperature  with  frequency  for  the  composite  problem 
for  ratios  of  the  diffusivity:  D2/D1  =0.1, 0.5,  5.0, 10.0. 
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Eigenvalue, 


frequency,  o)* 


(b)  Crack  located  in  the  glass 

Figure  24.  Stress  intensity  variation  with  frequency  for  composites  of  Glass  and  SiC. 


vn.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  conclusion,  we  have  shown  that  the  mathematical  techniques  presented  can  be 
used  to  evaluate  the  stress  intensity  factors  of  components  subjected  to  thermal  fatigue. 
The  applicability  of  the  evaluation  method  has  been  demonstrated  to  be  quite  broad. 
Conditions  of  asymmetric  temperature  distributions,  varying  geometries  and  differing  heat 
sources  can  be  simulated  through  the  use  of  this  method.  Because  of  the  wide  range  of 
applicable  conditions,  it  is  possible  that  we  can  use  these  results  in  laboratory  testing 
situations  to  approximate  the  contribution  of  the  thermal  fatigue  loading  in  thermal 
mechanical  fatigue.  Further,  we  can  extend  this  method  to  approximate  stress  intensities  in 
composite  materials.  An  example  is  the  case  of  periodically  occurring  cracks  in  fibers  of  a 
bimaterial  composite. 

The  results  presented  indicate  that  thermal  fatigue  loading  has  the  potential  to  cause 
significant  stress  intensities  in  general  cases  of  thermo-mechanical  fatigue,  and  the  method 
of  electric  resistance  heating  results  in  increasing  stress  intensity  factors  with  frequency. 
These  results  are  significant  enough  to  warrant  some  experimental  verification  of  these 
data.  A  much  needed  set  of  data  is  the  temperature  values  (magnitude  and  phase)  at  near¬ 
crack  and  far-field  locations,  induced  by  both  the  electric  resistance  and  hot  lamp  heating 
methods.  With  these  values  one  could  evaluate  the  magnitude  of  the  stress  intensities 
induced  in  actual  applications. 


97 


REFERENCES 


1 .  G.  C.  Sih,  "On  the  Singular  Character  of  Thermal  Stresses  near  a  Crack  Tip,"  Journal 
of  Applied  Mechanics,  29,  587-589  (1962). 


2.  A.  L.  Florence  and  J.  N.  Goodier,  "Thermal  Stresses  due  to  Disturbance  of  Uniform 
Heat  Flow  by  an  Insulated  Ovaloid  Hole,"  Journal  of  Applied  Mechanics,  27,  635- 
639  (1960). 


3.  G.  R.  Irwin,  "Fracture,"  in  Handbuch  der  Physik,  volume  VI,  S.  Fluegge,  ed.. 
Springer- Verlag,  551-589  (1958). 


4.  M.  L.  Williams,  "Stress  Singularities  Resulting  from  Various  Boundary  Conditions  in 
Angular  Comers  of  Plates  in  Extension,"  Journal  of  Applied  Mechanics,  19,  526-528 
(1952). 


5 .  Y.  Konishi  and  A.  Atsumi,  "The  Linear  Thermoelastic  Problem  of  Uniform  Heat  Flow 
Disturbed  by  a  Two-Dimensional  Crack  in  a  Strip,"  International  Journal  of 
Engineering  Science,  11, 1-7  (1973). 


6.  H.  Sekine,  "Thermal  Stress  Singularities  at  Tips  of  a  Crack  in  a  Semi-infinite  Medium 
under  Unifor  Heat  Flow,"  Engineering  Fracture  Mechanics,  7, 713-729  (1975). 


7 .  H.  Sekine,  "Thermal  Stresses  near  Tips  of  an  Insulated  Line  Oack  in  a  Semi-infinite 
Medium  under  Uniform  Heat  Flow,"  Engineering  Fracture  Mechanics,  9,  499-507 
(1977). 


8.  H.  Sekine,  "Crack  Problem  for  a  Semi-Infinite  Solid  with  Heated  Bounding  Surface", 
Journal  of  Applied  Mechanics,  44,  SZl-WL  (1977). 


9.  R.  Shail,  "Some  Steady  State  Thermoelastic  Stress  Distributions  in  the  Vicinity  of  an 
External  Crack  in  an  Infinite  Solid,"  International  Journal  of  Engineering  Science,  6, 
685-694  (1968). 


10.  B.  R.  Das,  "Some  Axially  Symmetric  Thermal  Stress  Distributions  in  Elastic  Solids 
Containing  Cracks  - 1.  An  External  Crack  in  an  Infinite  Solid,"  International  Journal  of 
Engineering  Science,  9, 469-478  (1971). 


98 


99 


11.  B.  R.  Das,  "A  Note  on  Thermal  Stresses  in  a  Long  Circular  Cylinder  Containing  a 
Penny-Shaped  Crack,"  International  Journal  of  Engineering  Science,  7,  667-676 
(1969). 


12.  K.  Herrmann  and  K.  Kuemmerling,  "A  Crack-Thermal  Stress  Problem  in  a  Doubly- 
Connected  Solid,"  Archives  of  Mechanics,  28, 171-188  (1976). 


13.  G.  B.  Sinclair,  M.  Okajima  and  J.  H.  Griffin,  "Path  Independent  Integrals  for 
Computing  Stress  Intensity  Factors  at  Sharp  Notches  in  Elastic  Plates,"  International 
Journal  for  Numerical  Methods  in  Engineering.,  20, 999-1008  (1984). 


14.  M.  E.  Gurtin,  "On  a  Path -Independent  Integral  for  Thermoelasticity,"  International 
Journal  of  Fracture,  15,  R169-R170  (1979). 


15.  W.  K.  Wilson  and  I.-W.  Yu,  "The  Use  of  the  J-Integral  in  Thermal  Crack  Problems," 
International  Journal  of  Fracture,  15, 377-387  (1979). 


16.  S.  Aoki,  K.  Kishimoto  and  M.  Sakata,  "Energy  Release  Rate  in  Elastic-Plastic  Fracttire 
Problems,"  Journal  of  Applied  Mechanics,  48,  825-829  (1981). 


17.  S.  Aoki,  K.  Kishimoto  and  M.  Sakata,  "On  the  Path  Independent  Integral-J," 
Engineering  Fracture  Mechanics,  13, 841-850  (1980). 


18.  S.  Aoki,  K.  Kishimoto  and  M.  Sakata,  "Elastic-Plastic  Analysis  of  Crack  in  Thermally 
Loaded  Structures,"  Engineering  Fracture  Mechanics,  16, 405-413  (1982). 


19.  A.-Y.  Kuo  and  P.  C.  Riccadella,  "Path-Independent  Line  Integrals  for  Steady-State, 
Two-Dimensional  Thermoelasticity,"  International  Journal  of  Fracture,  35,  71-79 
(1987). 


20.  Y.  C.  Fung,  Foundations  of  Solid  Mechanics,  Prentice-Hall,  Inc.,  p.  389  (1965). 


21.  D.  A.  Wilson  and  J.  R.  Warren,  "Thermal  Mechanical  Crack  Growth  Rate  of  a  High 
Strength  Nickel  Base  Alloy,"  ASME  Transactions,  paper  no.  85-GT-12,  30th 
International  Gas  Turbine  Conference  and  Exhibit,  Houston,  TX,  March  18-21, 1985. 


22.  F.  B.  Hildebrand,  Advanced  Calculus  for  Applications,  2nd  Edition,  Prentice-Hall, 
Inc.,  484-490  (1976). 


23.  R.  W.  Hombeck,  Numerical  Methods,  Prentice-Hall,  Inc.,  148  (1975). 


100 


24.  W.  H.  Press,  B.  P,  Flannery,  S.  A.  Teukolsky  and  W.  T.  Vetterling,  Numerical 
Recipes,  Cambridge  University  Press,  121-126  (1986). 


25.  W.  H.  Beyer,  ed.,  CRC  Standard  Mathematical  Tables,  27th  Edition,  CRC  Press, 
Inc.,459  (1984). 


26.  H.  Tada,  P.  C.  Paris  and  G.  R.  Irwin,  The  Stress  Analysis  of  Cracks  Handbook,  Del 
Research  Corporation,  Hellertown  PA,  p.  2.2  (1973). 


27.  N.  Marchand,  D.  M.  Parks  and  R.  M.  Pelloux,  "Ki-solutions  for  Single  Edge  Notch 
Specimens  under  Fixed  End  Displacements,"  International  Journal  of  Fracture,  31,  53- 
65  (1986). 


28.  P.  D.  Hilton  and  G.  C.  Sih,  "A  Laminate  Composite  with  a  Crack  Normal  to  the 
Interfaces,"  International  Journal  for  Solids  and  Structures,  7, 913-930  (1971). 


29.  T.  S.  (Zook  and  F.  Erdogan,  "Stresses  in  Bonded  Materials  with  a  Oack  Perpendicular 
to  the  Interface,"  International  Journal  of  Engineering  Science,  10, 677-697  (1972). 


30.  G.  D.  Gupta,  "A  Layered  Composite  with  a  Broken  Interface,"  International  Journal 
for  Solids  and  Structures,  9,  1 141-1154  (1973). 


31.  M.  Okajima,  Analysis  of  Tensile  Testing  Configurations  for  Assessing  the  Strength  of 
Butt  Joints,  Ph.  D  Dissertation,  Camegie-Mellon  University  (1985). 


32.  J.  P.  Dempsey  and  G.  B,  Sinclair,"  On  the  Singular  Behavior  at  the  Vertex  of  a 
Bimaterial  VJedge,”  Journal  of  Elasticity,  11, 317-327  (1981). 


33.  A.  R.  Zak  and  M.  L.  Williams,  "Crack  Point  Stress  Singularities  in  a  Two-Material 
^edge,"  International  Journal  of  Applied  Mechanics,  30, 142-143  (1963). 


34.  L.  G.  Loitsyanskii,  Mechanics  of  Liquids  and  Gases,  Pergamon  Press,  291-294 
(1966). 


35.  S.  E.  Cunningham,  J.  H,  Griffin  and  G.  B.  Sinclair,  "On  Stress  Intensities  Induced  by 
Direct  Resistance  Heating,"  International  Journal  of  Fracture,  33, 135-144  (1987). 


APPENDIX  A 


HEAT  SOURCE  INDUCED  BY  AN  ELECTRIC  CURRENT 
IN  A  CRACKED  PLATE 

Here  we  present  the  derivation  of  a  heat  source  in  a  cracked  plate  induced  by  electric 
potential  heating.  The  derivation  follows  the  method  of  Loitsyanskii  [34]  who  developed 
the  solution  to  the  analogous  problem  of  irrotational  flow  of  an  incompressible  fluid  with 
no  circulation  flowing  past  an  infinite  cascade  of  plates  (see  Figure  Al).  The  solution  to 
the  complex  velocity  of  the  fluid  is  derived  in  [34]  and  the  main  points  are  reiterated  here. 
Upon  integration  of  the  velocity  we  acquire  the  complex  potential  of  the  fluid.  The  electric 
potential  in  the  cracked  plate  is  analogous  to  the  real  part  of  the  complex  potential.  Once  the 
real  pan  of  the  complex  potential  is  found,  the  heat  source  is  obtained  by  taking  the  square 
of  the  gradient  of  the  electric  potential. 

We  consider  the  complex  coordinate,  2=x+iy.  According  to  the  geometry  of  Figure 
Al,  the  edges  of  the  plates  which  impede  fluid  flow  are  located  at  positions  in  the  complex 
plane: 


z  =  a  ±  2kW.  z  =  -a  ±  2kW,  k  =  0,1,2,... 


(Al) 


We  know  the  complex  velocity  of  flow  past  a  single  plate  in  an  infinite  medium  is 


V(z)  =  lu  -  i  v« 


z 


(A2) 


where  u*  and  v*  are  the  far-field  components  of  the  applied  velocity  in  the  x  and  y 
directions,  respectively. 

If  there  are  2k  +  1  plates,  then  z  in  the  numerator  of  (A2)  is  relaced  by  the  product 


Yl  (z  +  2kW)  =  z  J^(z2  -  4k^W^) 

k=-K  k=l 


(A3) 
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Replacing  similar  expressions  for  (z+  a )  and  (z-  a )  in  the  denominator,  taking  the  limit  as 
and  considering  the  identity 


k'hju 

k=i  ^  ^ 


=  sin 


Ttz 

2W 


(A4) 


gives  the  resulting  complex  velocity  of  fluid  flowing  around  an  infinite  cascade  of  plates: 


V(z)  =  u«  -  i  v^ 


Icos^-^ —  cos^-^^ 
V  2W  2W 


(A5) 


The  complex  potential,  0  (z) ,  is  obtained  by  integrating  the  complex  velocity: 


<I>(z)  =  Jv(z)  dz  =  u^z  -  ^  v,,^  log 


+c 


(A6) 


We  take  the  constant  of  integration  to  be  zero  without  loss  of  generality. 

The  complex  potential  can  be  written  as  the  sum  of  real  and  imaginary  parts  which 
represent  the  velocity  potential  and  stream  function,  respectively: 

<I>  (z)  =  <|»  (x,y)  +  i  V  (x.y)  ..j. 


In  the  electric  potential  analog  to  the  flow  problem,  the  velocity  potential,  <j»,  is  the  electric 
potential  and  I  =  -V  (|)/p  (where  p  is  the  electrical  resistivity)  is  the  current.  The  electric 
potential,  then,  is  obtained  by  finding  the  real  part  of  Equation  (A6)  when  the  far-field 
velocity  components  are  replaced  by  the  potential  drop,  <|>o.  across  the  plate  length.  If  the 
current  flow  is  assumed  normal  to  the  crack  this  gives  Uoo  =  0  and  Voo  =  <|>(/2L. 


♦  (x.y)  =  ^ 


(A8) 


Obtaining  the  real  part  of  the  above  expression  requires  some  algebra.  We  write 
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Z  =  cos  - CQS^-^^  =  Rc* 

2W  V  2W  2W 


e 


(A9) 


Then 


log  (Z)  =  In  R  +  i  0 


(AlO) 


To  find  R,  we  invoke  the  trigonometric  identity, 

cos||.  =  cos^cosh|^-ism^sinh|^  =  a-i3 


where  a  and  P  are  introduced  for  brevity.  In  addidon,  we  introduce  Oo: 


=  ao 


(All) 


(A12) 


Then 


Z  =  a 


ip  +  7(a-iP)^-o^ 


(A13) 


Letting 

^e‘’’=  (a-i  p)^  -oj 


(A14) 


then  the  square  rxxJt  portion  of  (A  1 3)  is  equal  to  ..for 


4  =  -J  (a^  -  P^  -  +  4  P^  *  "H  -  arctan 


2ap 


^2  o2  2 

a‘-p  -05 


(A15) 


Combining  these  expressions  gives 
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Z  =  a-ip±v^  e^'^  =  Re‘® 


(A16) 


Since  R2  is  equivalent  to  the  sum  of  the  squares  of  the  real  and  imaginary  parts  of  Z,  In  (R) 
is  equal  to 


In  (R)  =  1  In  {  [  Re  (Z)f  +  [  Im  (Z)  ]  ^ 


where 


Re(Z)  =  a±Jl  cosy 


and 


Im(Z)  =  -P±y^  sin^ 

Thus,  the  electric  potential  solution  is  given  by: 


(A17) 


(A18) 


(A19) 


0  =  ^-^  ln[  +  +  %±2yfl  (acosy  -psin-y)] 

2L  L  2  2  J 

The  plus  sign  in  ±  gives  solutions  for  positive  y  (+)  while  the  minus  sign  results  in 
solutions  for  negative  y  (-),  the  solution  being  symmetric  about  the  y-axis. 

We  obtain  the  heat  source  term  in  the  heat  conduction  equation  by  taking  the  square 
of  the  gradient  of  the  electric  potential  and  multiplying  by  the  electrical  conductivity  and 
dividing  by  thermal  conductivity: 


h  (x.y  )  =  y  (V  <|)  •  V  <J> ) 


a 

k 

. 

dyj 

(A21) 
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We  express  the  partial  derivatives  in  closed  form  and  compute  the  heat  source  with  the  aid 
of  a  computer  for  specific  points.  Near  the  crack  tip,  the  heat  source  behaves 
assymptodcally  as  1/r  where  r  is  the  distance  from  the  tip  of  the  crack.  Consistent  with  the 
condition  that  the  complex  velocity  approaches  a  constant  velocity  in  the  far-Held,  the  heat 
source  also  approaches  a  constant  value  at  far-field  points. 


Figure  Al,  Flow  past  an  infinite  cascade  of  plates 


APPENDIX  B 


GRIDS  FOR  FINTIE  ELEMENT  ANALYSIS 

Here  we  present  the  grids  used  in  the  finite  element  analysis.  The  number  of 
elements  and  degrees  of  freedom  for  the  homogeneous  grids  are  listed  in  Table  Bl. 
Figures  Bl  through  B3  illustrate  the  grids  used  in  the  homogeneous  problems.  Table  B2 
lists  the  number  of  elements  and  degrees  of  freedom  for  the  grids  used  in  the  composite 
analysis  and  the  grids  are  illustrated  in  Figures  B4  through  B6. 
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Table  B1 

Finite  element  grids  for  the  homogeneous  problems. 


Grid 

Number  of  elements 

Degrees  of  fircedom 

coarse 

72 

182 

medium 

288 

750 

fine 


720 


1562 
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Figure  Bl.  Finite  element  grid  for  homogeneous  stress  analysis,  coarse. 


Figure  B2.  Finite  element  grid  for  homogeneous  stress  analysis,  medium. 


■■■■■■■■■■■■■■■■■■■■■■■I 
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Figure  B3.  Finite  element  grid  for  homogeneous  stress  analysis,  fine. 
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Table  B2 

Finite  element  grids  for  the  composite  problems. 


Grid 

Number  of  elements 

Degrees  of  freedom 

coarse 

64 

162 

medium 

256 

578 

fine 

352 

794 
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Figure  B5.  Finite  element  grid  for  composite  analysis,  medium. 


Figure  B6.  Hnite  element  grid  for  composite  analysis,  fine. 


APPENDIX  C 


COMPLEMENTARY  EIGENFUNCTIONS  FOR  A  CRACK. 

Here  we  present  the  complementary  eigenfunctions  for  the  crack  in  an  elastic  plane. 
The  eigenfunctions,  established  in  polar  co-ordinates  with  the  origin  at  the  crack  tip,  are 
taken  fix)m  Sinclair  et  al.  [13]  and  are  specialized  to  the  crack  geometry.  The  symmetric 
forms  are: 


r 

■/2rt(l+K) 


cos^e 


Q*  - - If - --3/2  _1  cos -I-  6 --I- COS 

“  /Sfd+K)  L  2  2  2  2j 


'""ys  (1+K) 


TItc  (1+k) 

U  =_— I - 

-/27C  (H-K) 


/ 

[(i+K)cos|e-|«>s|] 

.[(i-ic)5in|e+|sin|] 


The  antisymmetric  complementary  fields  are  given  by; 


■flU  (1+k) 

isin^-lsin-^l 

2  2  2  rj 

1 

II 

\ 

[-^sin^  +-isin-|6] 

®®  /In  (1+K) 

[2  2  2  2  J 

X*  _  ^  r  ^ 

-^cos-®-  +-|-cos-|^ 

/2n  (1+K) 

2  2  2  2  . 
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IL  - - i - 

/27C  (1+k) 

u!=-— - 

®  yin  (1+K) 


where 


K  =  3  -  4v 
in  plane  strain,  and 


K  = 


3-v 

1+v 


in  plane  stress. 


isin|..(i+K)sin|e] 

■^COs|-  +  (j-K)  COs|a 


APPENDIX  D 


WEIGHTS  AND  ABSCISSAS  FOR  THE  GAUSS-TYPE  QUADRATURE 


1 

i 

0  ''  i=l 

Xi 

Wi 

n=l 

1 

0.333333333333 

2.000000000000 

n=2 

1 

0.115587109997 

1.304290309725 

2 

0.741555747146 

0.695709690275 

n=3 

1 

0.056939115967 

0.935827869145 

2 

0.437197852751 

0.721523146096 

3 

0.869499394918 

0.342648984758 

n=4 

1 

0.033648268067 

0.725367566757 

2 

0.276184313872 

0.627413291756 

3 

0.634677476235 

0.444762068907 

4 

0.922156608492 

0.202457072581 

n=5 

1 

0.022163568807 

0.591048449430 

2 

0.187831567652 

0.538533438620 

3 

0.461597361496 

0.438172725032 

4 

0.748334628387 

0.298902698301 

5 

0.948493926288 

0.133342688617 

n=6 

1 

0.015683406607 

0.498294091627 

2 

0.135300011655 

0.466985073077 

3 

0.344942379427 

0.406334853446 

4 

0.592750127732 

0.320156657087 

5 

0.817428013267 

0.213878651991 

6 

0.963461278703 

0.094350672773 

n=8 

1 

0.009027377026 

0.378901220910 

2 

0.079300559811 

0.365206830090 

3 

0.209779368616 

0.338313038790 

4 

0.381771053397 

0.299191977633 

5 

0.570635820162 

0.249257942511 

6 

0.749317378547 

0.190317023365 

7 

0.892221974214 

0.124507047877 

8 

0.978914210162 

0.054304918824 
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n=10 


n=12 


n=16 


n=24 


1 

0.005856308437 

0.305506774261 

2 

0.051886393980 

0.298345972945 

3 

0.139656240743 

0.284192218637 

4 

0.260985093682 

0.263377276898 

5 

0.404564284766 

0.236389063923 

6 

0.557011314600 

0.203860239634 

7 

0.704117292400 

0.166553483153 

8 

0.832171652087 

0.125344096668 

9 

0.929241876580 

0.081202859601 

10 

0.986304414519 

0.035228014278 

1 

0.004103285523 

0.255876390693 

2 

0.036526421504 

0.251674912694 

3 

0.099251890030 

0.243340945856 

4 

0.188176807259 

0.231011336107 

5 

0.297484581452 

0.214888540232 

6 

0.420025381681 

0.195237304208 

7 

0.547783818960 

0.172380323064 

8 

0.672403257001 

0.146692962822 

9 

0.785732486514 

0.118597169831 

10 

0.880359134936 

0.088554877635 

11 

0.950095757826 

0.057062777258 

12 

0.990397602845 

0.024682459600 

1 

0.002333630564 

0.193080177029 

2 

0.020872147684 

0.191277440159 

3 

0.057258441734 

0.187688798162 

4 

0.110136769181 

0.182347757392 

5 

0.177536897897 

0.175304186009 

6 

0.256947517675 

0.166623848454 

7 

0.345409811309 

0.156387791574 

8 

0.439627699909 

0.144691588218 

9 

0.536090655003 

0.131644445553 

10 

0.631204502056 

0.117368186957 

11 

0.721425343258 

0.101996118525 

12 

0.803391614274 

0.085671796045 

13 

0.874049370825 

0.068547725826 

14 

0.930766209806 

0.050784130618 

15 

0.971430051691 

0.032548789462 

16 

0.994535210151 

0.014037220019 

1 

0.001048475472 

0.129475393625 

2 

0.009409911669 

0.128932328872 

3 

0.025992648096 

0.127848477169 

4 

0.050518761473 

0.126228384573 

5 

0.082577199139 

0.124078846320 

6 

0.121630668223 

0.121408878332 

7 

0.167024640562 

0.118229679397 

8 

0.217998322465 

0.114554584201 

9 

0.273697405472 

0.110399007400 

10 

0.333188384403 

0.105780378970 
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11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

n=48  1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 


0.395474202728 

0.459510963069 

0.524225422743 

0.588532981170 

0.651355857686 

0.711641155150 

0.768378506728 

0.820617010337 

0.867481167625 

0.908185562375 

0.942048039521 

0.968501196288 

0.987102203204 

0.997543524928 

0.000264932423 

0.002382707518 

0.006609282769 

0.012926746332 

0.021308325370 

0.031718499509 

0.044113151374 

0.058439753552 

0.074637591201 

0.092638019352 

0.112364753820 

0.133734194491 

0.156655779610 

0.181032369576 

0.206760658608 

0.233731612547 

0.261830930934 

0.290939531402 

0.320934054339 

0.351687385674 

0.383069195575 

0.414946490776 

0.447184178188 

0.479645637414 

0.512193299734 

0.544689231106 

0.576995716721 

0.608975844627 

0.640494085946 

0.671416869237 

0.701613146562 

0.730954948851 

0.759317928231 

0.786581884998 

0.812631277024 

0.837355709413 

0.860650402358 


0.100718071108 

0.095233316985 

0.089349121713 

0.083090165887 

0.076482702132 

0.069554445130 

0.062334455666 

0.054853019417 

0.047141521679 

0.039232320915 

0.031158631446 

0.022954469158 

0.014655107803 

0.006306692105 

0.065101228985 

0.065032237428 

0.064894327428 

0.064687645137 

0.064412409588 

0.064068912464 

0.063657517789 

0.063178661542 

0.062632851194 

0.062020665173 

0.061342752247 

0.060599830842 

0.059792688273 

0.058922179916 

0.057989228301 

0.056994822130 

0.055940015234 

0.054825925452 

0.053653733451 

0.052424681471 

0.051140072011 

0.049801266445 

0.048409683585 

0.046966798172 

0.045474139317 

0.043933288878 

0.042345879784 

0.040713594309 

0.039038162280 

0.037321359255 

0.035565004632 

0.033770959729 

0.031941125805 

0.030077442054 

0.02818188354S 

0.026256459134 

0.024303209342 
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38 

0.882416635207 

0.022324204200 

39 

0.902562164861 

0.020321541070 

40 

0.921001616744 

0.018297342462 

41 

0.937656846722 

0.016253753851 

42 

0.952457272490 

0.014192941582 

43 

0.965340173199 

0.012117091008 

44 

0.976250956557 

0.010028405486 

45 

0.985143394070 

0.007929108677 

46 

0.991979831560 

0.005821463636 

47 

0.996731426993 

0.003707921578 

48 

0.999379104174 

0.001593584131 

i 


